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We present a new pseudospectral code, bamps, for numerical relativity written with the evolution 
of collapsing gravitational waves in mind. We employ the first order generalized harmonic gauge 
formulation. The relevant theory is reviewed and the numerical method is critically examined and 
specialized for the task at hand. In particular we investigate formulation parameters, gauge and 
constraint preserving boundary conditions well-suited to non-vanishing gauge source functions. Dif¬ 
ferent types of axisymmetric twist-free moment of time symmetry gravitational wave initial data are 
discussed. A treatment of the axisymmetric apparent horizon condition is presented with careful 
attention to regularity on axis. Our apparent horizon finder is then evaluated in a number of test 
cases. Moving on to evolutions, we investigate modifications to the generalized harmonic gauge 
constraint damping scheme to improve conservation in the strong field regime. We demonstrate 
strong-scaling of our pseudospectral penalty code. We employ the Cartoon method to efficiently 
evolve axisymmetric data in our 3 + 1 dimensional code. We perform test evolutions of Schwarzschild 
perturbed by gravitational waves and by gauge pulses, both to demonstrate the use of our blackhole 
excision scheme and for comparison with earlier results. Finally numerical evolutions of supercrit¬ 
ical Brill waves are presented to demonstrate durability of the excision scheme for the dynamical 
formation of a blackhole. 


I. INTRODUCTION 


This is the first in a series of papers about the nu¬ 
merical treatment of collapsing gravitational waves using 
a new pseudospectral code developed for the problem. 
In the early 1990s critical phenomena were discovered 
in gravitational collapse [I], in spherical symmetry, with 
general relativity minimally coupled to a massless scalar 
field. One aspect of the finding was that, amazingly, 
the critical solution dividing the formation of a blackhole 
from dissipation of the field, was unique, in the sense 
that if one takes any one parameter family of initial data, 
with the parameter controlling somehow the strength of 
the data, and tunes this parameter to the threshold of 
blackhole formation, one finds that the same solution is 
always obtained, regardless of the family! Shortly there¬ 
after similar phenomenology was reported in axisymmet¬ 
ric, vacuum general relativity [2], or in other words in the 
collapse of gravitational waves. Since then multiple stud¬ 
ies have been performed to reproduce this finding, albeit 
with different initial data and numerical approaches, but 
without success. Perhaps most strikingly, in [3], numeri¬ 
cal evidence of a different critical solution was presented. 
Even if one completely accepts the available evidence for 
criticality in vacuum collapse, this obviously begs the 
question whether or not the naive expectation of unique¬ 
ness of the critical solution in axisymmetric, rather than 
spherical, collapse holds. 

Roughly speaking there are two types of code being 
used used in 3d numerical relativity. The first uses 
the moving puncture method El E]> which consists, in 
essence, of a clever choice of evolved variables and gauge 
conditions, normally treated numerically by finite differ¬ 
encing. Secondly is the pseudospectral method, most 
prevalently used with a first order generalized harmonic 
formulation of general relativity by the SpEC code [6]. 
Recently, we presented a study of the collapse of grav¬ 


itational waves using the moving puncture method |7], 
in part to establish how close to the critical regime one 
can get with this standard approach. The conclusion be¬ 
ing; not very. Major difficulties included the formation 
of coordinate singularities and a lack of accuracy. There¬ 
fore one would like to tackle the problem using the pseu¬ 
dospectral approach to establish what can be achieved 
in that setting. We have thus developed a new pseu¬ 
dospectral code along the lines of SpEC, specializing the 
continuum and numerical method as much as possible 
towards the problem of vacuum gravitational collapse. 
The present paper represents the first outcome of this 
maneuver. Herein we describe the formulation of GR em¬ 
ployed, our boundary conditions, the code, calibration of 
the method, our initial data, our approach to axisyrn- 
metric apparent horizons, plus a suite of validation tests 
for gauge waves, gravitational waves, blackhole and col¬ 
lapse spacetimes. Throughout we compare our results 
carefully with those in the literature. We aim to give a 
body of evidence for the correctness of the method that 
the reader will find compelling. With this out of the way, 
in subsequent papers we turn to the problem of critical 
collapse. A summary follows before the main text. 

In section [II] we look at a slightly modified version of 
the first order generalized harmonic formulation of [5]. 
We consider constraint preserving, radiation controlling 
boundary conditions, paying special attention to the con¬ 
straint preserving boundaries. By considering the reflec¬ 
tion of outgoing waves in the linear approximation we 
ultimately suggest modified conditions that should re¬ 
duce spurious reflections caused by the use of constraint 
damping. We also suggest alternative gauge boundary 
conditions. 


Next, in section III we outline the bamps code, includ¬ 
ing our carefully constructed cubed-sphere grids, which 
avoid clustering of grid-points in unfortunate positions 
of the domain. For the discretization we employ a pure 
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Chebyscliev approach. We also discuss our ‘octant’ sym¬ 
metry implementation, the crucial patching-penalty ap¬ 
proach for communicating data between neighboring co¬ 
ordinate patches, and finally the boundary implementa¬ 
tion. In the follow-up section[lV]we complete the presen¬ 
tation of the penalty method by computing the penalty 
parameters appropriate for the semi-discrete system. 

Given the difficulties in the literature in reproducing 
the results of [ 2 ] it seems necessary to solve the problem 
in axisymmetry before moving to examine the collapse of 
fully 3d waves without symmetry. In our moving punc¬ 
ture gauge study [7] a major disadvantage in using the 
BAM code was that 3d grids were employed to evolve ax- 
isymmetric data. In section[V]we present our approach to 
evolving axisymmetric spacetimes with the bamps code, 
for which we employ the Cartoon method [9] to reduce 
from the standard bamps 3d domains to a plane, by using 
the Killing vector to evaluate any angular derivatives. We 
discuss various flavors of axisymmetric moment-of-time- 
symmetry initial data and their numerical construction. 
These initial data sets are evolved in a forthcoming study. 
We also give a detailed description of our formulation of 
the apparent horizon conditions in axisymmetry. To the 
best of our knowledge this is the first time that the reg¬ 
ularity conditions on the symmetry axis have been care¬ 
fully taken care of. This is important in later work as the 
search for apparent horizons will be our key diagnostic 
tool. 

The next three sections ( VT|| V III ) contain a write-up of 
our development and validation tests. The tests include 
evolutions with the proposed gauge boundary conditions, 
which we find are helpful when using large gauge source 
parameters, as desired. They also include runs compar¬ 
ing the fully 3d, octant symmetry and Cartoon evolu¬ 
tions, demonstrating that the various symmetry setups 
are well-behaved. In the evolution of single blackholes 
we test different gauges and boundary conditions, and 
following [Sj, look at evolutions in which the blackhole is 
perturbed by a gravitational wave injected through the 
outer boundary. Our results are in good agreement with 
the earlier studies. We then examine the evolution of su¬ 
percritical Brill waves, where, after the formation of an 
apparent horizon the run is continued after interpolation 
onto an excision grid, as used to evolve a single black- 
hole, which is needed to evolve data with a horizon for 


long-times. Finally we conclude in section IX 


II. THE GENERALIZED HARMONIC 
FORMULATION AND BOUNDARY 
CONDITIONS 

A. GHG, constraints, boundary conditions 

The evolution system: We use the first order reduc¬ 
tion of the generalized harmonic formulation with several 
free parameters. The full reduction from the second or¬ 
der Einstein equations is presented in detail elsewhere 0 


so here we give only a brief overview to establish our no¬ 
tation. Throughout the paper in continuum equations 
we use the latin a,b,c... for four dimensional indices, 
but for spatial indices, with the exception of n 

and s, whose meaning when used as indices will be de¬ 
scribed shortly. Greek indices are used to refer to the po¬ 
sition in a state-vector, grid indices, or where otherwise 
needed. We start from the vacuum generalized harmonic 
formulation in second order form, 

Rab = V (o C b) + 7 4 r c a 6 C c - \^9ab9 cd T e cd C & 

- 7o \n( a C b ) - g a bn c C c ] , (1) 

for the unknown spacetime metric g a b with Christof- 
fels T c ob . The constraints of the system are C a = 
g bc T abc + H a = 0, plus the standard Hamiltonian and 
momentum constraints of GR. The gauge source func¬ 
tions H a are freely specifiable, provided that they do not 
include derivatives of the metric, which would affect the 
principal part of the PDE. The terms involving 70 are in¬ 
cluded so as to damp away high-frequency constraint vi¬ 
olations m- The parameters 74 and 75 control whether 
or not the constraint addition made in the construction 
of the formulation is done either with the covariant or 
the partial derivative, or some combination. The latter 
choice has the effect of simplifying the constraint sub¬ 
system. In the code we use a first order reduction by 
introducing the variables d> iab and n ab . The equations of 
motion are, 

fttVab ft dig a b ^Rat> T 7l ft ^iab 7 
fttft&iab ft^ ftj^iab 7y0iTl ab T 72Q((A ab -p ^CXTl 71 ^icdX^ab 
+ a^ k n c ^ij C ^ ka b 7 

dtRab = ft l diR ab - a^di^jab + lil2ft l C iab 

+ ‘2‘ag cd (j l i<& ica <&j db — n ca n db - g ef T ace T bd f ) 

- 2 a(V {a H b) + 74 T c ab C c - i 75 g ab r c C c ) 

- |an c n d n cd n ab - an c ^T\ ci ^ 3 ab 

-p « 7 o [2(5 (a^G) gab^l ] C c , ( 2 ) 

with shorthands to be defined momentarily. The formu¬ 
lation here agrees with that of [ 8 ] except for the inclusion 
of the 74 and 75 parameters. We will either take the new 
parameters to vanish, or choose 74 =75 = 1/2. The 
lapse and shift are denoted a and ft 1 respectively. The 
unit normal to the spatial slices of constant coordinate 
time t is written n a . When the normal is contracted with 
a tensor we sometimes use the abbreviation S an = S ab n b , 
and likewise for the arbitrary unit spatial vector s a . The 
induced metric on the slice is written 7 ^. In matrix no¬ 
tation this system can be written as 

d t u» = A k » v d k u v + S'*, (3) 

with u M = (g ab 7 n ab , <l>i ab ) r , and principal matrix, 

/(l + h)ft k 0 0 \ 

= 7 i72/3 fc ft k , (4) 

\ 72cu5f —aS k ft k ) 
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and containing all source terms. We use the short¬ 
hand for the Christoffel symbols under the first order 
reduction, 

h abc — r y\b\^i\c)a 2*7 a^ibc A 'bl(^b^-c)a 2^ a B-6c 5 (b) 

and will frequently use the abbreviation T a = g bc T a bc . 
The system is symmetric hyperbolic, having the same 
principal part as a particular first order reduction of the 
wave equation. The characteristic variables are given by, 

u° ab = dab , 

U ab = n af> T sl ®iab ~ 72 9 ab , 

U Aab = tfA'&iab , (6) 

with the projection operator q J t = 8 J i — s :l S{. and speeds, 

u S = (1+ 7 i)/ 3 s , v ± =/3 s ±a , v? = (3 s , (7) 

respectively. For future reference let us also note that a 
convenient way to transform to the characteristic vari¬ 
ables is to write u a = T~ la p , where here the indices 
represent the position in the state-vector u a and where 
the similarity matrix is, 

( i o o\ 

T~ 1A a= 72 1 f , (8) 

M -72 Is* w 

V 0 0 

which has left inverse But note however 

that T~ la p 7^ 6 a j}. The strength of this represen¬ 
tation in practical terms is in avoiding special cases in 
the numerical implementation, like for example s x = 0, 
in the characteristic decomposition. 

Gauge source functions: For the gauge source func¬ 
tions H a we choose, 

H a = Vl fog (~—^ n a - ^§7 aiP 1 - (9) 

\ a J or 

Our convention differs from that of both El and 3j in 
a trivial normalization of the spatial part with respect 
to the lapse function. Writing the resulting gauge condi¬ 
tions in terms of the lapse and shift we get, 

d t a = - a 2 I< A rj L a 2 log + fPdia , 

dt? = a 2 (3) F* -ad 1 a- nsF A /3 j d j p i , (10) 

with K the trace of the extrinsic curvature and (3) r l the 
contracted Christoffel symbol of the spatial metric. Be¬ 
fore blackhole formation for the scalar functions 774,, 77 s 
we choose, 

7 1 l = f) L a q , 77s = fj s a r , (11) 

with ?7L,77s,g,r some constants. By default we 
choose p = 1 and q = r = 0, which naturally maintains 


the shift damping term even if the lapse function is close 
to zero, in contrast to the standard condition employed 
in SpEC El. which takes r = 1 . Since we wish to study 
near-singular gravitational effects in the computational 
domain and avoid run-away growth of the shift vector 
this seems reasonable. We will report in later work on 
adjustments to these choices when evolving near-critical 
data. When evolving blackholes by excision we follow HD 
taking instead r = 1, although so far we have not found 
it necessary to use the log 2 form of 77 l- 

The constraint subsystem: The first order reduced 
harmonic constraints are, 

Ca = Ha A 7 VQija - + n b U ab - \n a g bc Ii bc . 

( 12 ) 

The terms without H a are simply r a . In these vari¬ 
ables the vacuum ADM Hamiltonian and momentum 
constraints can be expressed as 

2 G nn = 7«7 fc, (0*$ tfI - d k $ Hj A r a jk T ail 

- T a ijT ak i ) , 

7 i G na 7 ^ (cl[j T 2 dj*hkin 2 di*&jkn 

2 ^j[^j] nn + 7 *&mk\j*&i]ln 
+ 2T an[l T a k]j ). (13) 

As stated above we use a subscript n to denote contrac¬ 
tion with the normal vector 7i a , but with the convention 
that di stands for the partial derivative, but with any 
such contraction outside of the derivative. We can put 
the Hamiltonian and momentum constraints together as 
a four-vector of constraints, 

M a = G an . (14) 

Working with the first order system creates the reduction 
and closely related ordering constraints, 

Giab — dig ab 4r } i a b — 0 i 

G’ijab — di*hj ab dj&i ab — 2d^jC{^ a b 9 • (lb) 
The constraints C a and Ci ab evolve according to, 

d t c a = (1 A 71 )FdiC a - 71 P l diC a +aG a 

A (74 - "/5)an a T b C b - a(2q 4 - l)r b an C b 
+ 2q 0 an b n( a C b) A 0 : 7^7 kl $i kn Cij n n a 
-a 1 i a C ijn [\g bc & bc + & nn ] 

- 7i72/3* ( \g cd Ci cd n a - C ina ) , 
dtCiab = ( djC iab A 71 diCjab) A a [(1 + 71 )d l gj„ C J ab 

72 Giab A ab Gij n A ^Gi nn Tl ab j , (19) 

where we have introduced the constraint, 

G a = 2 M a A (n a7 ib - 7 \n b mc b - T c lb C c ) 

+ 72(S c a 'f ib -y >c Y a )C ibc , (17) 
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and where the notation di means take the partial deriva¬ 
tive, and afterwards replace all first derivatives of the 
metric with the reduction variable &i a b- Up to lower 
derivatives in the contraints we find, 

d t G a « PdiGa + on^didjCa - cn ik 'y li d l c ijka 

+ \a^al ll g cd diC locd , (18) 

where « denotes equality up to non-principal terms, the 
remainder having been suppressed for brevity. The equa¬ 
tion of motion for Cij ab is readily derived by taking 
derivatives of that of Ci a b . Notice that the parameter 72 
serves to damp the reduction constraint. In the descrip¬ 
tion of 0 the equivalent reduction variable is called F a , 
with, including 74 and 75 in the natural way, 


as 

c m on = J d 3 x^y (s ab F a F b + S ab C a C b + Y j 6 ab C ia C jb 

+ 7 ij S ac S bd C iab C jcd + 7 « 7 kl 6 ac 8 bd C ikab C jlcd ). 

(23) 

The gravitational wave degrees of freedom: In vacuum 
the Weyl scalars T 0 , T 4 can be expressed as, 

*o = m A m B [L^ bd AB l a l c Rabcd ], 

*4 = m A m B [± (F)6d 4 B k a k c R abcd ] , (24) 

respectively. Here we have introduced the null tetrad 


Fa = G a - (1 - 74 ) (n„ T b - 2 T b an )C b - 75 n a Y b C b , 

(19) 


in our variables. The difference is not substantial, being 
only that G a appears slightly more naturally in the sec¬ 
ond order form of the equations. Note that in (19), the 
final term contains a piece which is simply the Harmonic 
constraint in the pure harmonic case, but will act as a 
non-zero coefficient otherwise. 

First order reduction of the constraint subsystem: 
Following H], a first order reduction of the constraint 
subsystem is formally introduced by defining the new 
variable Ci a with, 


C ia = 7 3 k dj®ika ~ §7; ° a g cd dj^icd + diH an - \n a g cd dill cd 

+ diH a + + h jk ®F C ®iknn a 

- 7 j W m ^jia^ikm + ^ icd ll be n a (g cb g de + \g be n c n d ) 

- $ icn Il ba (g bc + \n b n c ) + 57 2{n a g cd - 2 5 c a n d )C icd . 

( 20 ) 


The principal part of this formal reduction is given by, 

d t c a « 0 , 

dtG a * FdiG a +Off*diC ja , 
dtCiaKFdjCia + adiGa, 
dtCiabttil + riFdjCiab, 

dtCijab ~ /3 k d k C ijab . (21) 

The characteristic variables of the constraint subsystem 
are then found to be 

4 = F a T C sa , C° = C a , 

C Aa = d A.Gia , c[ab = ^iab , 

C ijab = Gijab , (22) 

with speeds /3 s =F a, 0, /3 s , (1 + 7 i)/ 3 s and /3 s respectively, 
where we use upper case latin indices to denote those 
projected by q a b . A suitable norm of the constraint vio¬ 
lation is given by the constraint monitor which is defined 


r = -G(n a +s a ), k a = ^(n a -s a ), 

m a = Uj ( v a + iw a ) , m“ = ^ (v a - iw a ), (25) 

with s a ,v a and w a mutually orthogonal unit spatial vec¬ 
tors, and the projection operator, 

I (P)cd c. 1 „cd „ 

-L ab — Q (aQ b ) 2 ^ ^ a b 

= m^ a m b ^fh^ c m d ' > + fh^ a m b yn^ c m d ^ . (26) 

In terms of the first order GHG variables we can express 
the principal part of the Riemann tensor as, 

ddabcd ~ 7'^ off *Yj b [c^fd] 4 b@i^ja[cTd] T Cl r , dj 11 fr[c7d] 

7Z^5jH a [ c 7ci] + 7 a^i^d b [c^d] 7 b^i^d-a[c^ / ct\ 

- n a 7 u di(j) jb[c n d] + n b Y J di(j) ja [ c n d] 

- 7i72 n a n k d k g b [ c n d ] + 'yi"f 2 n b n k d k g a [cn d ] 

- 727 l adig b [ c n d ] + 727 \dig a [ c n d ] ■ (27) 

Of course this expression is unique only up to constraint 
additions. Note that upon contraction with _L (P ) and l 
to form the Weyl scalar To, and after a single addition 
of Cij ab , we naturally form a projection of the incoming 
characteristic variable d s u A b . This is used in the con¬ 
struction of the boundary condition. The spatial vec¬ 
tor s l is taken to be the unit spatial normal to the bound¬ 
ary. 

Boundary conditions: At the outer boundary we need 
to control incoming constraint violation, gauge pertur¬ 
bations and physical radiation. By default we initially 
impose, 


F a +C sa + lC a = 0, (28) 

on the constraint subsystem assuming that the charac¬ 
teristic variable cf is always incoming. These conditions 
are essentially those of [ 8 j, with just the additional 1 /r 
term. Other conditions for this variable will be motivated 
and tested in what follows. The remaining constraint 
subsystem characteristic variables may or may not be in¬ 
coming, and are dealt with on this basis as described in 
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section III D[ but always according to the same prescrip¬ 
tion. For the gravitational wave degrees of freedom we 
choose, 


the effect of the different constraint preserving boundary 
conditions let us consider now the subsystem without the 
reduction. We have, 


^0 = 90 ) 


(29) 


the lowest order member of a cascade of conditions on 
incoming radiation nana, with given data <? 0 - Examin¬ 
ing (241 it is obvious that this is equivalent to setting, 

^ P)M AB(n C Rated) = 1 (P) % , (30) 


which is in practice how the conditions are implemented. 
For the remaining gauge degrees of freedom we choose 
either the improved gauge boundary conditions of El, 


_L 


ft )Cd d t [Kd + (72 -r l )g cd ] = 0, 


(31) 


or the alternative, 


X 


(G)cd 

ab 


dsV'cd d)\ H - ^2^sed 

+ 7,-1 ( u td - 2n (c H d) + i2g C d) 


= 0 , 


(32) 


\ 7 b Y ba = —R ab C b , 

Y ba = V b C a + 2 l4 T c ab C c - (74 - 7 5 )3a b r c C' c 

- 27 0 n ia C b ) . (35) 

The shorthand Y ab and the variable G a that follows will 
be related to quantities present in the first order reduc¬ 
tion of the GHG formulation shortly. We can equivalently 
express this as, 

n b d b C a = G a - (274 - l)T c ab n b C c - (74 - "i b )n a T c C c 
+ 270 n b ri( a C b ), 

n b d b G a = 7 bc V b [V c C a + 2 l4 T d ac C d - (74 - i b )ga C T d C d 

- 2 lo n (a C c) \ + (n b V b n c )[V c C Q + 2 l4 T d ac C d 

- (74 - l b )gaeV d C d - 270 n (a C' c )] + T c ab n b C c 

+ RatC b , (36) 

where the variable, 


(G) 

with given data q y cd ', which we will often take to van¬ 
ish, and where the overbar derivative notation has the 
same meaning as in equation ( |17[ ). These conditions are 
similar to the ‘freezing’ gauge boundary conditions em¬ 
ployed in [ 8 ] , but taking into consideration the discussion 
of gauge reflections given in [T4], and constructed so that 
the conditions are naturally applied to metric compo¬ 
nents (in ADM form) and their derivatives, but excluding 
contributions from the gauge sources. We will typically 
try to choose the given data to be fixed in time, and 
such that initially the time derivatives vanish for these 
quantities. Here we have introduced the gauge projection 
operator, 

-Lib )cd = l(ak b) l {c k d ) + k a k b ri d - 2k {a q b) W . (33) 

The above boundary conditions are implemented in 
bamps using the Bjprhus method [15] as in SpEC. Details 
of the method are explained in section |IIID| For com¬ 
pleteness here the constraint projection operator X (C ) = 
/- ±( p ) - ±( G ) is, 


G a = n b Y ba = 2M a + (n al ib - ^ a n b )V t C b , (37) 


is used to allow for the most convenient form of these 


expressions, and the final term of (36) is in fact of sec¬ 


ond polynomial order in the constraints because of the 
vacuum field equations ([!]). Different choices of the con¬ 
straint addition parameters 74,75 result in different be¬ 
havior in terms of growth of the constraint fields. It 
is also obvious that different choices of these parameters 
can simplify the constraint subsystem, the natural choice 
apparently being 74 = 75 = 1 / 2 . 

Linearization: Let us linearize and consider the be¬ 
havior of a set of fields that satisfies these equations on 
a fixed constraint satisfying background. We start with 
equation (35) and use the tetrad consisting of the null 


vectors l a ,k a 1 m a ,m a defined in (25) to decompose the 


first index of Y ba . From this we obtain, 


V 6 k b l c Y ca + 4 k c Y ca - m b m c Y ca - m b m c Y ca = 0, 


(38) 


Y { a ? Cd = I Qab q cd - 2 l (a q b) ^ + l a l b k c k d , (34) 

and also plays an important role in the imple¬ 
mentation of the boundary conditions, as they are 
again naturally written in the form _L^ cd d s uf d = 
transverse derivatives. 


B. Constraint preserving boundary conditions and 
damping 

Generalized harmonic constraint subsystem: We al¬ 
ready saw the constraint subsystem of the first order re¬ 
duction of the GHG system. But to get a better idea of 


for the linearization, where we are free to use the no¬ 
tation C a for the linearized violation because the con¬ 
straints are satisfied in the background. 

Boundary conditions: Taking the standard setup at 
the outer boundary so that s a , used in the construction 
of the tetrad, denotes the outward pointing spatial unit 
vector normal to the boundary. Restricting our atten¬ 
tion to boundary conditions that contain at most one 
derivative of the constraints, geometrically the most nat¬ 
ural choice seems to be l b Y ba = 0 . In the first order GHG 
language these conditions are, 

Ga + S7 s C a + 2 74 r c as C c 
+ (74 - 75 )s a r & C b - 70 n a c s = 0. 


( 39 ) 
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Whereas, discarding the first order reduction, those 
of (28) are instead, 


Ga + v s c a + T c as C c - (2 74 - l)T c an C c 
- (74 - l 5 )n a r b C b + \C a = 0. (40) 


With either conditions one might guess that the 
choice 74 = 75 = 1/2 reduces reflections from the bound¬ 
ary, especially when using a non-harmonic r a = — H a 7 ^ 0 
gauge. Incidentally this choice also makes the two condi¬ 
tions almost coincident. Suppose all derivatives of C a , G a 
tangent to the boundary vanish, and that the background 
is flat. Then we can analyze the solutions in a plane wave 
approximation. 

Mode solutions on flat-space: When linearized 
around flat-space this system takes the form, 


the remainder is rather of order 0(70 C a w _1 ). There is 
some freedom in expressing these conditions in the first 
order GHG language, but we choose, 

Ga + VsCa + 2 74 r C as C' c + (74 — 7 5 )s a r b Cfc 
+ 2 7o7a fc Cfc ~ 7o n a ( C n + ^C s ) + ^ C a = 0 . (45) 

The conditions ( |39[ ) can be similarly rewritten. A similar 
analysis can be performed using the pure gauge subsys¬ 
tem presented in m, but we currently find that existing 
gauge boundary conditions are sufficient for our needs, 
so we do not present these calculations here. Tests with 
the various boundaries are presented in section [VT| 


III. THE BAMPS CODE 


UC a - 270 d b (n( a C b )) = 0 . 

The right-travelling mode solutions are, 

_ ,, .sft+iaix 1 s sft+iuix 

On — Pl e 1 + P 2 e , 

/~ii _ i sft+iux 

O — P2 e 1 

with eigenfrequencies, 


(41) 

Having discussed the continuum system in the previ¬ 
ous section, we now discuss details of our numerical im¬ 
plementation of the GHG system. For this we present 
the bamps code, which uses a pseudospectral method on 
cubed-sphere grids. The basic idea of the code is based 

(42) on SpEC B], but in many details, such as the actual 
grid implementation and the outer boundary treatment, 
bamps differences are present. 


= -|7o - 
st = -7o - i 


4co 2 - 7 q 


oA- 7o A 


(43) 


A very desirable property for our boundary conditions 
would be that they absorb outward going waves per¬ 
fectly, that is, without reflection. With this motivation 
high-order derivative boundary conditions on the gravita¬ 
tional wave degrees of freedom have been studied mm, 
and implemented in the SpEC code |16 in order to ab¬ 
sorb higher spherical harmonics of the Weyl scalar 4/ 4 . 
In the current context absorption means that outgoing 
mode solutions, those associated with an s + , lie in the 
kernel of the boundary conditions. This is only the case 
if we switch off the damping 70 = 0. Since the low order 
spherical harmonics might be expected to dominate in 
the gauge and constraint subsystems, optimizing against 
this phenomena may be more important than using high- 
order conditions for the gauge and constraint subsystems 
whilst neglecting the damping terms. 

Remainder of mode solutions: Substituting these 
mode solutions into the boundary conditions (40), or the 
natural geometric conditions 


each after appropriate 
linearization, and expansion at large frequency u> gives 
remainders of order 0 ( 706 * 0 ), indicating that neither is 
the optimal that can be obtained by adding source terms 
to the constraint boundary conditions. Taking instead, 


(dt + d s + 70 )C n + \"ioG x — 0 , 

ipt + d 3 + \"io)Ci = 0, (44) 


A. Grid setup 

Grid types: The numerical domain on which we solve 
the evolution equations in bamps is either a cubed-ball 
or a cubed-shell grid. Each type is built up of multiple 
deformed cubes. Each patch is described by two fun¬ 
damental overlapping charts. In local coordinates cc, y 
and 5 it is a rectangular box x [—1,1] x [—1,1]. 

In global Cartesian coordinates x , y and z the cubes are 
transformed and rotated in such a way that when added 
together they build the desired domain. We give a de¬ 
tailed description in the following. The cubed-ball-grid 
includes the origin and has a spherical outer boundary. 
It consists of 13 coordinate patches: 

The central cube: is centered around the origin and 
ranges from —r cu to r cu in the global Cartesian co¬ 
ordinate directions. 

The transition shell: transfers the grid from the inner 
cube grid to a spherical shell with radius r cs . It 
contains six patches. 

The outer shell: consists of six patches which extends 
the grid with additional cubed-shells up to the 
outer grid boundary at r ss . 

The cubed-shell-grid is an excision grid, meaning that 
it does not include the origin. It is a special case of 
the cubed-ball-grid, consisting only of the six outer shell 
coordinate patches. 
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Cubed-sphere coordinate transformation: The coordi¬ 
nate transformation used in bamps to construct the grids 
introduced above relies on the so called “cubed sphere” 
construction. It was introduced in [IS] and first ap¬ 
plied in the context of numerical relativity in num. 
Since then this idea was implemented in multi-patch ap¬ 
proaches (!?TH24| . In contrast to many of the earlier ex¬ 
amples, the numerical method of bamps does not require 
overlapping grids, which simplifies the discussion. In m, 
the coordinates are constructed by considering great arcs 
parametrized by equidistant angles. Such angle coordi¬ 
nates are used in hhm], while [221 Hi] use an interme¬ 
diate set of coordinates also given in m that does not 
have the equidistant angle property. In bamps the latter 
type of coordinates are employed. The concrete coor¬ 
dinate transformation is the following. First, the local 
coordinates of each patch are transformed to temporary 
global coordinates 

x x_ x_ . . 

Xt = Vt= -V , Zt = (46) 

s s s 

This patch, which is orientated in positive x direction, 
will later be referred to as the master patch. From here, 
cyclic permutation is used to rotate the patches to their 
location in the sphere. The denominator s depends on 
where the coordinate transformation happens. For the 
patches of the outer shell it is 

s=(l + f + z 2 ) 1/2 , (47) 


we choose the number of subpatches to be J\f cn x A f cn . In 
Fig. a we show a 2d sketch of the bamps cubed-ball grid 
subdivided into subpatches. 

Discussion: It is straightforward to specify a mapping 
between a rectangular master patch and a cubed sphere, 
although some book keeping for the different patches and 
different types of shell transitions is involved. It may be 
useful to examine different such mappings in terms of a 
numerical quality criterion, say the size of the Jacobian, 
and to minimize the distortions associated with the co¬ 
ordinate transformation. 


B. Numerical method 

Spatial discretization: bamps uses the method of lines 
with a standard ODE integrator to integrate in time. The 
right-hand-sides are approximated using a pseudospec- 
tral method. We use a linear transformation to map 
the local coordinates of each subpatch x* into a unit 
cube x* = ( x,y,z) T £ [— 1, l] 3 . We discretize the sub¬ 
patch by choosing Gauss-Lobatto collocation points in 
each dimension, for example, 

7 T 

x a = - cos ( N _ 1 a) , (50) 


In the transition shell its definition includes a transition 
function A 


s(A) = (r 


1+2A \!/2 


+ A (y 2 + Z 2 )) 


A = 


Xt - Xn 


(48) 


This coordinate transformation is constructed to transi¬ 
tion from the inner cube to the outer shells. Note that 
this transformation is uniform along the 3d diagonals, 
where the distance between inner and outer shell bound¬ 
ary is smallest. This significantly improves the time¬ 
stepping restriction in the transition shell. 

Subpatches: Each coordinate patch can be further di¬ 
vided into subpatches. Subpatches are helpful for in¬ 
creasing resolution, and form the backbone of the par¬ 
allelism of bamps. Each master patch can be split 
into J\f x x M v x M z subpatches with coordinates 


with a = 0 , • • • , N x — 1, and similarly in the other direc¬ 
tions. The number of of grid points N depends on the 
patch location in the grid. The central cube is discretized 
with N cn x N cu x N cu points. The radial directions of 
the transition and outer shell are filled with N cs and N ss 
points respectively. The number of angular points we 
chose to be the same as in the central cube to assure 
that we have matching grids. In Fig. [I] we show on the 
right the Gauss-Lobatto discretization of a subpatch. 

Basis expansion: On the collocation points we ex¬ 
pand all evolution fields u in each dimension in a spectral 
basis using Chebyshev polynomials T ra (x) 


N x — 1 

u(x a , yp , Zs') ^ [ c n fyp, z§)T n (x a ), (51) 

71=0 


x l £ [xg + k 1 Ax 1 , Xq + ( k l + l)Ax*], (49) 

with Ax* = { Xl ^° , jj^) and fc* = 0 ,...,Mi - 1. In 
practice we ensure subdivisions are made in such a way 
that subgrids of two neighboring patches match, and that 
neighboring patches and subpatches share grid-point po¬ 
sitions on their respective boundaries. This is necessary 
because our current penalty-communication method does 
not deal with interpolating penalties. Concretely we split 
the inner cube into M cu x A/" cu x M cu subpatches. The 
transition and outer shell are divided in M cs or M ss sub¬ 
patches in the radial direction. For the angular direction 


and analogously in the remaining two directions. We use 
the pseudospectral approach and store not the expansion 
coefficients c x ,c v ,c z but the function values u a p$ at the 
collocation points x l a p S . 

Derivatives: The spatial derivatives of the evolution 
fields are computed by a matrix multiplication. For ex¬ 
ample in the x-direction we have 

N x -1 

(dxu) a j35 = 'y [ D arL U n /3S ( 62 ) 

71=0 







Central box 



FIG. 1: The left part of the diagram gives a two dimensional sketch of the bamps cubedball grid layout. The ball is built up of 
several transformed cubes. These patches can further be divided in subpatches. In the example shown we have J\f cu = 3,Af cs = 2 
and Af sa = 1. On the right is shown that each subpatch is covered by Gauss-Lobatto grids ranging from —1 to 1 in local 
coordinates. 


with the Gauss-Lobatto derivative matrix, 




2(JV x -1) 2 + 1 
6 

go (-ir +g 

qp x a —xp 
-xp 

2(i-£|) 
2(jVx-l) 2 + l 
6 


a = p = 0 

/3 

a = P = !,-■■ ,N X — 2 

a = p = N x - 1 


(53) 


where q a = 2 at boundary points and q a = 1 elsewhere. 
In practice we do not compute diagonal terms of the 
derivative matrix by the analytic formulas stated above 
but use the identity 


Ate —1 

D aa = - Y, £, »- ( 54 ) 

n=0,n^a 


This negative-sum-trick maps a constant function explic¬ 
itly to zero and is known to give the derivative matrix 
better stability as regards rounding errors 1251 . In prelim¬ 
inary experiments we found that this gives slightly more 
accurate derivatives, but have not studied the influence 
on the accuracy of the simulations presented later in the 
paper. 

Filtering: We find that a crucial ingredient for numer¬ 
ical stability is the use of a filter against high-frequency 
growth. For this we follow EH 3 exactly. After every full 
time-step we apply the filter in each dimension. The fil¬ 
ter is easily implemented as a matrix multiplication. For 
example, in the x direction, we filter the function values 

by, 


n 


N x -1 

E 


t^nl 37 : 


(55) 


with the filter matrix 

•Aa/3 = 51 S an e- 36(n/n — r A n0 , (56) 

n 

where n max = N x — 1 and S a p and A a) 3 are the Chebyshev 
synthesis and analysis matrices respectively. 

Time integration: We integrate the fields forward in 
time using a 4th order Runge-Kutta scheme. Unless oth¬ 
erwise stated we fix the time-step, At = ^Ar m j n , with 
Ax m i n being the minimal Cartesian spatial grid spacing 
of the whole domain. Empirically we find that this choice 
for the time step always leads to stable numerical evolu¬ 
tions, in the sense that increasing resolution results in 
smaller errors. We have not not carried out a stability 
analysis of the fully discrete system. 

BAMPS Octant grid: When evolving octant sym¬ 
metric data in bamps, it is possible to only evolve one 
eighth of the cubed ball grid. This saves computational 
and memory costs. In the bamps octant mode we choose 
an odd number of subpatches A 4 u and a odd number 
of grid points N cu and reduce the numerical domain 
to x > 0, y > 0 and z > 0. This means that all sub¬ 
patches containing one of the Cartesian axes are cut in 
half along them. In these patches we use the symmetry 
conditions to construct special matrices which compute 
the derivatives and Liters. 


C. Patching boundary conditions 

To glue all subpatches together, we have to impose ap¬ 
propriate conditions at the connecting boundaries of the 
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subpatches. For this we apply the penalty method as 
is described in [26ll28] . The main idea of this method 
is to add penalty terms for each incoming characteristic 
variable at the boundary to the right hand side of the evo¬ 
lution equations. We use the characteristic variables of 
the evolution system to formulate boundary conditions. 
On the boundary surface we define the outward point¬ 
ing spatial normal vector s l . The characteristic variables 
of the evolution system are given in equation ([6]) with 
speeds ([7|. In vector notation we write 

U ab \ 

4 • ( 57 ) 

U AabJ 

Incoming characteristic variables to the subpatch bound¬ 
ary have positive speeds. On these we want to impose the 
condition that they are equal to the outgoing character¬ 
istic variables of the neighboring patch. Table |T] summa¬ 
rizes all incoming and outgoing characteristic depending 
on the lapse function a and the shift in s l direction, (3 s . 
As an example, let us now consider the boundary be¬ 
tween two patches, patch L and patch R , and the case 
—a < (3 s < 0. With respect to the spatial normal vec¬ 
tor s‘ at the boundary pointing outwards of subpatch L 
and inwards in subpatch R , the incoming characteristic 
variables of L are the outgoing ones of R. In the chosen 
case, u+ h are incoming to L and outgoing of R. We want 
to impose the condition, 

< b L =< b R . (58) 

Multiplying the first order GHG evolution equations from 
the left with the matrix of eigenvectors T -1 a p, we obtain 
evolution equations for the characteristic variables. 



D. Outer boundary implementation 


At the spherical outer boundary of the domain we use 
the Bjprhus method 2} H5j to impose the constraint, 
physical and gauge conditions given in section [TTJ As for 
the patching boundaries, we impose conditions on the in¬ 
coming characteristic to the boundary surface. However, 
this time instead of adding penalty terms we modify the 
right hand side of the evolution equations at the bound¬ 
ary in such a way that the boundary conditions are satis¬ 
fied. We define the outward pointing spatial normal unit 
vector s l and use the projection operator q 3 , t = S 3 i — s 3 Si 
to split the principal part of the evolution equation in a 
part normal and tangential to the boundary surface 

d t u^ A k ^(s k s 3 + q^dX 

= A s ^d s u" + A A » v q%d B u 1 '. (61) 


Expressed in characteristic variables the normal part is 


dtu & ~ T - 1 * li A" i v T' $ T- l &td a ue = A a& $ d a J (62) 


The matrix A sa ^ is a diagonal matrix containing the 
characteristic speeds. At the outer boundary we assume 
that the absolute value of the shift (3 s is always smaller 
than the size of the lapse a. This leads to two cases to 
be considered. 

Case —a < /3 s < 0: In this case the incoming char¬ 
acteristic at the outer boundary condition is u + . Ac¬ 
cording to section [TT] we impose the following boundary 
conditions, which we give here only schematically: 


1. One of the constraint preserving boundary condi¬ 
tions (p8|), (|39|) or (451), 


J_ (G) d s u++ P (C) + NP^ =0. (63) 


d t u &L =T~ 1& fl A k ^ u d k u vL +T~ 1& /1 S^ (59) 

Here the d again denotes that the similarity ma¬ 
trix T~ la ^ stands outside the partial time derivative. 
At the boundary we now add a penalty to the right hand 
side of the evolution equation of the incoming character¬ 
istic. This is often called weakly imposing the boundary 
condition, 


d t uj b L = T - 1 +„A k ^d k u v ab L + T- x +^ 

+ p(< b R ~<b L ) ■ ( 60 ) 


Afterwards we use the inverse transformation to get back 
to the evolution equations enhanced with the necessary 
penalty terms at the boundary. These are also the equa¬ 
tions we implement in the code. We treat all six bound¬ 
aries of the subpatches independently from each other. 
This means that on the edges we have to consider penalty 
contributions from two and on the corner from three di¬ 
rections. The size of the penalty parameter p can be 
derived from an energy estimate of the semi-discrete evo¬ 


lution system. This we present in a section IV 


2. One of the gauge boundary conditions (31) or 
which become either, 

_L (g) d t u+ + P (G) + iVP (G) = 0, 

_L (G) d s u+ + P (G) + iVP (G) = g (G) . 


(64) 


3. The physical boundary condition (29), 

-L (P) d s uj b + P (p) + ATp( p ) = 9 ( p ) . 


(65) 


Here we labeled principal terms with derivatives tan¬ 
gent to the boundary P lx ) an( t non-principal terms with 
NPC). At the boundary surface we project the evolu¬ 
tion equation of the incoming characteristic u + into the 
constraint, the physical and gauge part. 

dtuj b « v+(±^ cd + ±4 cd + cd )d s uj d . (66) 

All three parts have to be replaced using the boundary 
conditions. We do this by subtracting the conditions 
from the bulk right hand side D t , 

d t u+ b = D t u+ b - u + (Conditions) a& , (67) 
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/3 s > a > 0 

a> /3 s > 0 

p a = 0 

-a < /3 s < 0 

/3 s < -a <0 

'U'ab 

0 

zero 

zero 

zero 

zero 

zero 

Kb 

f3 s -a 

incoming 

outgoing 

outgoing 

outgoing 

outgoing 

“aft 

/3 s + a 

incoming 

incoming 

incoming 

incoming 

outgoing 

U Aab 

r 

incoming 

incoming 

zero 

outgoing 

outgoing 


TABLE I: Incoming and outgoing characteristic variables to a subpatch boundary with spatial normal vector s 1 depending on 
the gauge variables. 


with the special case (31) treated in the obvious way. 
Transforming back this modified right hand side leads to 
modified evolution equations at the boundary. 

Case 0 < /3 s < a: In this case also the characteris¬ 


tic u Aab is incoming. As described in |8], we impose the 
additional constraint preserving boundary condition 


ds u A ab ~ q B A d B ®sbc = 0 , ( 68 ) 


by subtracting it from the evolution equation of u Aab 


IV. ENERGY ESTIMATE FOR PENALTY 
FACTOR 

In this section we derive an estimate for the right 
choice of penalty factor at the patching boundaries of the 
bamps domains. The actual technical implementation of 
the patching condition was already described in subsec¬ 
tion |III C[ The following calculation is based on the one 
presented in 128 . However we present it for a general hy¬ 
perbolic system in curvilinear coordinates, albeit under 
rather restrictive assumptions. 


d t' u Aab = D t^Aab ~ (Condition Aab ). (69) 


A. The continuum case 


After we have modified the right hand sides at the bound¬ 
ary we transform back to the evolution equations for the 
primitive fields. 


E. Code implementation details 


Code structure: The bamps code is written in the C 
programming language in a modular fashion. The code 
is designed in such a way that the technical layer is sep¬ 
arated from projects for solving physics problems. Inside 
physics projects we use a Mathematica script, MathToC, 
which translates equations written in tensor notation into 
C code. As a standalone program we have developed an 
axisymmmetric apparent horizon finder, AHloc, which is 
typically used to search apparent horizons in bamps gen¬ 
erated data at the post-processing step. It is also possi¬ 
ble to run the finder in a daemon-like mode in which it 
searches horizons in data of a running instance of bamps. 
We describe the apparent horizon in subsection VE| 
Parallelization: bamps is programmed to run in paral¬ 
lel on several computing nodes using the message passing 
interface (MPI). The N s u b subpatches of a bamps grid 
are distributed on M MPI processes as evenly as pos¬ 
sible. This means that each process has to handle at 
least n = subpatches. As in general the total num¬ 

ber of grids is not divisible by the number of MPI pro¬ 
cesses without remainder, IV S ub mod M processes have 
to take care of one additional grid. In practice we choose 
the number of MPI processes in such a way that the num¬ 
ber of processes which have to compute one grid less is 
minimized. 


We view the GHG system as a general symmetric hy¬ 
perbolic system of partial differential equations, but sup¬ 
press all non-principal terms, and work in the linear, con¬ 
stant coefficient approximation, so we have, 


dtu* = AP^d p u\ 
where, in matrix notation, 


p<Ex,y,z, 


9ab \ 

/(l + 7i)/? fc 

0 

n ab 1 

, A p K = 7l72 /?* 

pk 

*&iab J 

V 'rz a5 i 

-at i 


(70) 


-cry 

ok 


ik 


(71) 


For clarity we suppress the state vector indices p, v. For 
this system there is a symmetrizer H such that HA p s p 
is Hermitian for every unit spatial vector s p . The energy 
of the system is, 


E 2 


iv 


dV {u*Hu). 


(72) 


with the volume form dV = dx d y dz v /q. As discussed in 
section III A each subpatch of bamps has a set of global 
Cartesian coordinates x l = (x, y , z) and a set of local 


coordinates x l = ( x,y,z ). The Jacobian J~ = ^ trans¬ 
forms between the two charts. To formulate boundary 
conditions at the patching boundaries which control the 
energy in the patch, we study the time derivative of the 
energy, using the evolution equations we replace the time 
derivatives by spatial derivatives, 


d t E 2 


d Vd p [vtHA p u] . 


(73) 
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In the constant coefficient approximation we can com¬ 
mute the determinant of the three metric in the volume 
form with the partial derivative and end up a divergence 
in flat Cartesian coordinates 

d t E 2 = J dxdydzdp \dHA p Uydy] . (74) 

In the next step we change to the patch local coordinates 
x, y and z, 

d t E 2 = J dV [u'HAtuyfy det Jj\ 

= J dxdy d z dp& p . (75) 

Here we have defined := y/ldet J~ and the flux <f>P = 
id H APuy/^f. Now we integrate over all boundary surfaces 
of the patch, 

ii ii 

8 t E 2 = J Jdydz$% = „ 1 + J J dxdz&ll^! 
- 1-1 - 1-1 

1 1 

+ JJ dxdy<t> s \l = _ r (76) 

-l -l 

At a boundary surface, for example x = const, we can 
write the unit normal vector as, 

t = {y k djxd k x)~* 7 il dix = ld l x, (77) 

S_ v ,✓ 


and 2 + 1 split the spatial metric jij ■, 


7 ij = 




(78) 


The relationship between the determinant of 7 + and the 
metric in the boundary surface q is, a/t = lyfq- We 
rewrite, 

<fr x = 4 ddpx = d HA p uly/^dp x = HA a u , (79) 


and express the time derivative of the energy as the sum 
of boundary surfaces integrals over the fluxes >I> S , 

ii ii 

dtE2 = / / d ' 4 »> 4 ‘l»--.+// d - 4 « 4 ‘ll=-i 

- 1-1 - 1-1 

1 1 

+ JJdA^'\l_ 1 . (80) 

-l -l 

The area element is dA yz = y/q dy dS. The fluxes can 
be rewritten in terms of incoming and outgoing charac¬ 
teristic variables at the boundary surface. The system 


is symmetric hyperbolic. Therefore the principal symbol 
has a full set of Eigenvectors which we write as columns 
of the similarity matrix T a . With the inverse of this ma¬ 
trix, T s _1 , we transform the vector of evolution variables 
to the characteristic variables of the system v = T s _ 1 u. 
The flux expressed in the language of characteristic vari¬ 
ables is, 

-I s = u t (T s " 1 ) t TlHT s (: T s )~ 1 A s r, T~ 1 u = v^HAsV. 

jjt fj A„ v 

(81) 

The diagonal matrix A s contains all the speeds of the 
characteristic variables 


A s = 




(82) 


Where we have ordered the characteristic variables in 
such a way that we group all incoming with positive 
speeds A j and outgoing with negative speeds —A jj. In 
this partition it follows that 





and with this 


= v\Hih.iVi - vljHuAnVn . 


(83) 


(84) 


If all integrands in (801 are negative semi-definite, the 
energy of the system does not grow over time. For the 
boundary conditions we use the ansatz vj = kvji + 3 , 
which means that at the boundary surface we set the 
incoming characteristic variables equal to a linear combi¬ 
nation of the outgoing characteristic variables plus some 
given data g. Choosing the matrix idn small, we obtain, 


l> s = (. g f + u} J K t )^/A/(KU// + g) - vljHnAnVu 


< g'H r Arg + v' u 


id HjAjk — HjjAjj 


VII . 


(85) 


The first term only depends on the given data. As we 
are free to choose it we have full control over this term. 
The second we can make negative again by choosing d k 
sufficiently small. 


B. The semi-discrete case 


In this subsection we carry out the energy estimate for 
a semi-discrete system. In our case this means that we 
discretize the evolution variables in space using Gauss- 
Lobatto collocation points according to equation (50). 
The semi-discrete evolution equations are 


d t u a ps = A^[d p uU 5 = A p [JP] a/3 s[dpu\ a ps ■ ( 86 ) 
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The energy of this system is defined using Gauss-Lobatto 
quadrature with the appropriate integration weights 

w ee> w j8> W «J 

E 2 = • ( 87 ) 

o/38 

Again we compute the time energy of the system, 
with £j a ps = oj a u)poJsWl\oi 05 i using the inverse product 
rule to write, 


with the penalty matrix 

P « - (T o) • (94) 

and 8v a ps = [v BC ] a p 5 — \pf\ a ps, with v BC the desired 
boundary data. The time derivative of the energy splits 
into two parts 

d t E 2 =d t E^ + d t E 2 pen . (95) 


dtE — ^ '^ apS^p 
o(38 


J J o(38^ a fi8-^ > ^ J ot/38 


( 88 ) 


and transform to local coordinates. For this we assume 
that dp^/'yaps = 0 and obtain, 


dtE — 'y ' ^apsdp 
o(38 

Using an expansion in Legendre polynomials we can use 
the summation by parts property to write, 


UaasEaPS^UaPsVy, 


o(38 


(89) 


The first part is the contribution from the bulk, 

dt Abuik = ^ 2^05 b/]0/3(5 [H] 0/35 [Af ]ops [vi ]006 
05 

- y ^U 06 [v\i] 005 [Hllhps [Afj]o,95 [vil]o 05 ■ 

05 

(96) 

The second part changes the time derivative of the energy 
because of the additional penalty terms in the evolution 
equation at the boundary, 


dtE — y ' ups ‘^ j a p^E a 05 A~ u a 0s\fy a 


05 


T y ^ 03 a5 'U’apsEapS-A^Ufxpsy/^tf‘ 


a.(38 


08 


+ u lx 05 E a 0 sA U a 0 S\/^ a p S 

a(3 


N x -1 

CK=0 
Ny-1 
0=0 
N z -1 
5=0 


(90) 


^■^pen E ^O05([uq 05^ H 0 psT 0 p S Pps8v 0 p5 

05 

+ ^Vops]^pls T O05 H O05U O p5) ■ (97) 

By inserting the identity TT = I into the appropri¬ 
ate places we transform the state vector u to the vec¬ 
tor of characteristic variables. Then multiplying out the 
penalty matrix and rearranging leads to, 


As in the continuum case we introduce the normal out¬ 
ward pointing s l vector at the boundary and write, 


dtE — y ' UJps '^' a psE a psAJ*ctpS'U’apS 
05 

T y [ LJgS ^ a psE[ a psAP[Sp] a p5U a p§ 
08 

^ ^ ^a/3 ^ a /3$Ha(38AP[Sp\ a {38'U' 0i /38 
a(3 

with tips = \Zqcops- We define the flux 

*&a05 — 'U'apsEapd \Sp\a05'U’ap5 , 


N x -1 
a=0 
Ny~ 1 
0—0 
N z -1 
5=0 


(91) 


(92) 


and transform it to characteristic variables in the obvi¬ 
ous way. This gives us for the semi-discrete case the 
analogue expression for the time derivative of the energy 
at the boundary (80). In case of patching the boundaries 
between two subpatches we apply the penalty method to 
impose boundary conditions. For simplicity we restrict 
ourselves to the a = 0 boundary. For each incoming 
characteristic variable we add a penalty term to the right 
hand side of the evolution equations, 


dtU a ps — A \Jp\ot05 [^p^]a/35 T ^a, 0 [7s]/35-F/35 8 v a ps ■ 


(93) 


dtEp Cn = Pp5Coop5{[v BC ]lp S [H I ] 0 ps[v BC ]op5 
05 

- [ v Aops[Ei]o05[vi]o05 ~ [Ho/ 35 [Hi\o 06 [Ho/ 35 ) • 

(98) 

In total the change of energy at the boundary surface is, 

dtE 2 = y^IHoH^A/ ^ Ppswops)[Hi] 0 ps [u/] 0 ps 
05 

- E t A 3 < 5 b//] 0/35 [Ell] 005 [A//] 0/35 [vil] 0/35 
05 

+ '^2p05^O05[v BC }lp S [Hl}o05[v BC }o05 
05 

- y^P/35^0/35 [Hq/35[H]0/35[Ho/35 • (99) 

05 


We now consider two neighboring subpatches which we 
label L (for left) and R (for right). Let us assume they 
have a common boundary at a = N — 1 for the left patch 
and a = 0 for the right patch. For each subpatch we 
can write down the change of energy as in equation (991. 
As boundary conditions we set the incoming characteris¬ 
tic variables of one patch to be the outgoing one of the 
neighboring grid, 

v bc = V II ! 


— R 

V BC ~ V II > 


( 100 ) 
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and demand that the change of energy of the sub patches 
in time due to the patching boundary is not growing. 
Sufficient conditions for this are given by, 


P/38 


Wo /38 


Pps 


W(jV-l)/3<5 


( 101 ) 


In bamps we use these penalty parameters, but our dis¬ 
cretization is made with Chebyschev rather than Leg¬ 
endre polynomials, the equations we solve are not linear 
with constant coefficients and nor are the Jacobians map¬ 
ping from the master coordinates to our global Cartesian 
coordinates constant. Therefore it is to be determined 
empirically that the implemented method is in an appro¬ 
priate sense stable. These facts may contribute to the 
necessity of employing the filter ( 561 . 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 

x axis 

FIG. 2: The apparent horizons for centered Brill data 
with A = 11.82 and A = —5.3 and for pure plus polariza¬ 
tion data with A = 2.381 and A = —2.28. 


V. AXISYMMETRIC CONSIDERATIONS 

Although bamps is a fully 3d code we are often inter¬ 
ested in evolving axially symmetric data, which requires 
special attention for efficient treatment. In this section 
we collect together the relevant developments undertaken 
for axisymmetric initial data, apparent horizons and time 
evolution with the bamps code. 


oblate. In the initial data an apparent horizon can first 
be found at A = 11.82 with horizon mass M // = 4.8. 
For geometrically oblate data an apparent horizon can 
first be found at A = —5.30 with mass Mh = 4.4. The 
pseudospectral method we use to solve the constraints 
is discussed a little more in section |VD| Our apparent 


horizon search is explained in VE 


A. Brill wave initial data 


B. Pure plus polarization wave data 


Brill wave initial data is described in detail in many 
other sources. For completeness we give a bare-bones 
summary to highlight the particular choices that we 
make. 

Metric ansatz: Following 12131323 , we start from a spa¬ 
tial metric of the form, 

d^ 2 = ') i jAx' l dx^ = 4 / 4 [e 2l? (d p 2 + dz 2 ) + p 2 d(j) 2 ] , (102) 

in cylindrical polar coordinates, and take the extrinsic 
curvature to vanish. Note that the assumption of confor¬ 
mal flatness in the p-z sector of the metric can be made 
in axisymmetry without loss of generality. Under these 
assumptions the momentum constraints are trivially sat¬ 
isfied and the remaining Hamiltonian constraint takes the 
form 


D 2 T = - 


T 

J 




We then make the parametrized ansatz, 


(103) 


Metric ansatz: Observers distant from a compact ob¬ 
ject see gravitational waves in the form, 

d l 2 = dr 2 + r 2 ( 1 + h+)d0 2 + r 2 sin 2 61(1 - h+)d(/) 2 

+ 2r 2 sin# h x d6d(j ), (105) 


with the wave polarizations h+ and h x small perturba¬ 
tions of the Minkowski metric. This suggests modifying 
the ansatz (102) to 


dP = dr 2 + r 2 (e 2q dd 2 + e~ 2q sin 2 9d(f 2 ) 


(106) 


so that if we choose the seed function small and centered 
far from the origin we will have initial data that rep¬ 
resents a pure plus polarization gravitational wave. One 
could similarly make an ansatz for pure cross polarization 
waves, or indeed make other choices completely like [ 33 ] 
which we have also implemented and tested. 

The constraints: Again we start with moment of time 
symmetry initial data, so the remaining constraint takes 
the form, 


q{p, z) = Ap 2 e- [ ( p - po)2+[z - Zo)2] . (104) 

for the seed function q(p, z) and solve for T with bound¬ 
ary conditions T = 0 for asymptotic flatness at spatial 
infinity. This ansatz is the same as that studied in a num¬ 
ber of other studies GMaEDGE]. We call data with A > 0 
geometrically prolate, and that with A < 0 geometrically 


AV=^'SfR. (107) 

The conformal Ricci scalar is, 

R =~2 [ e ~ 2q - 1 ~ (r d r qj 2 ] - 1 3 dg(sin 3 dd g e ~ 2q ), 

r 2 r 2 sin 9 

(108) 
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and the Laplacian of the conformal metric is, 

1 e 2q 

A4' = -^d r (r 2 d r ^) + — ^-dg(e~ 2q sin 2 9dg^f) . (109) 
r 2 v 7 sin 9 v 7 

Once more we impose the obvious boundary conditions 
for asymptotic flatness at spatial infinity, and choose the 
seed function, 

q(r,9 ) = Ar 4 sin 2 9e~f r 2 - 2 rposine+p oJ, (110) 


which makes the conformal metric regular on axis. 

Apparent horizons: Taking centered data with A < 0 
we first find an apparent horizon at around A = —2.28, 
with mass Mh = 5.47. Looking for apparent hori¬ 
zons in centered data when A > 0, we find the curi¬ 
ous result that there is a region [2.381,2.568] in which 
apparent horizons are first found. Curiously, in the 
range [2.569,3.006], the data again seemed to be horizon¬ 
less. Continue at A = 3.007 we find horizons again up 
to A = 3.750 where we stopped our search. We searched 
for horizons using the resolution A A = 0.001. A closer 
look at the data at the boundaries of the ‘horizonless’ 
region shows that the shape of the horizon is very nearly 
not a ray-body, and we expect that our method simply 
can not find the horizons in this range of amplitudes (see 
section YE). We expect that this could be remedied by 


implementing an offset in p in the parametrization of the 
surface similar to that in z which we already have, but we 
leave this improvement for the future. The first apparent 
horizon for this data, found at A = 2.381, is plotted in 
Fig. [2] It has a mass of Mh = 4.8. 


C. Teukolsky wave initial data 


Initial data for numerical relativity: Teukolsky 
waves [34] [35] are an exact solution to GR linearized 
around flat-space, and were used as a seed function in jj], 
the first numerical study of the critical collapse of grav¬ 
itational waves, in the construction of full solutions to 
the constraints. In particular the waves were taken to be 
centered at some r 0 7 ^ 0 , with a radial width much less 
than ro, with an l = 2 , m = 0 spherical harmonic depen¬ 
dence, and mostly incoming. Since we are restricting to 
moment of time symmetry data, we can not satisfy the 
last of these conditions, but we expect that if the waves 
are placed at some sufficiently large ro then they will 
initially be weakly self-interacting, and roughly one half 
of the wave will simply propagate outwards. One could 
use the ansatz made in the Teukolsky wave initial data 
to construct incoming boundary data, but we leave this 
for future work. The construction of these data are well- 
described in [36] and were employed in [7]. See also m- 
Therefore here we want only to describe a subtlety that 
was overlooked in both of these references. 

Regularity of the conformal metric: Let us consider 
the ‘polar’ Teukolsky data. A similar discussion holds for 


axial data. The conformal metric for the solution of the 
Hamiltonian constraint is, in spherical polar coordinates, 


Irr = 1 + afrr , jrff = bf r g T , 

lee = (1 + cfee - a)r 2 , 

lu = (1 - cfgg + a f^)r 2 sin 2 9, (111) 

with the remaining components vanishing. Here we have 
already restricted the ansatz by removing terms that van¬ 
ish for l = 2 and m = 0 spherical harmonics. The angular 
functions f rr , f ge , f r e and /^ are, 

f rr = 2 — 3 sin 2 9 , f r e = — 3 sin 9 cos 9 , 

fee = 2 — 3sin 2 9 , /^ = 3sin 2 0 - 1, (112) 


whilst the remaining radial functions a, 6 , c are con¬ 
structed according to the recipe of [36], so that, 


a = 3 


b = — 


° 4 


F&> 3 FW 3 F 


r 2 

i?(4) 


r* 

3^(2) 

y3 

2F^ 


6 FW 

r 4 

9F( 2 ) 


6 F 

y*5 

21FW 


(113) 


21F 


In this expression we have used the shorthand, 


pi n ) 


d n F(x)~ 

dx n 


' d n F(x) ' 

dx n 


(114) 


and finally the seed function is F(x). In [7] the seed 
function was taken to be, 

F(x) = (e-^ x+ro)/a]2 + e -^ x - r °^^ 2 ) , ( 115 ) 

2 <t V / 

with p = 1. For local flatness however it is necessary |3T] 
that the combinations, 

cos 2 9 j rr + r~ 2 sin 2 9 jee — sin 29 ^ r g , 
r ^ 1 cos 9 "f rr — r~ 3 cos 9 ^gg + r ^ 2 sin ^ 1 f? cos 20 7 ^ , 
sin 2 9 'frr + r~ 2 cos 2 9 j eg + r~ 2 sin 2 9 7 ^ + r~ l sin 9 7 ge , 
r~ 2 Irr + r ~ 4 tan -2 9 7 gg - r~ 4 %</, + 2r ~ 3 tan -1 9 7 r g , 


of the metric components are regular functions of z = 
r cos 9 and p 2 = r 2 sin 2 9. Therefore one may worry about 
the high powers of r ^ 1 present in the recipe. This worry 
is justified, because for the particular seed function (115) 
the resulting conformal metric has a conical singularity 
at the origin, since for example the latter combination 
diverges like r~ 5 as r —> 0. Therefore the seed function 
used in [7] is not suitable to construct regular gravita¬ 
tional wave initial data. Regularity is obtained if one 
chooses instead takes p = 9. The fact that such a high 
power of x is required in the seed function illustrates the 
depth of the singularity that was present beforehand. 
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Comments on the numerical results o/PL with Teukol- 
sky initial data: Since the seed function (115) gives rise 
to an irregular conformal metric, the corresponding ini¬ 
tial data evolved in [7( were in principle wrong, as more 
careful convergence testing may have revealed. Should 
we then discard those results? Probably not; since the 
parameters taken were r 0 = 2 and a = 1/2, both the seed 
function and its derivatives were highly suppressed at the 
origin, and therefore in practical terms it is unlikely that 
the leading error in the simulations was caused by this 
problem. We will not attempt however to evolve the older 
data with bamps, not least because constructing irregu¬ 
lar data with our spectral elliptic solver is troublesome. 
These issues do not affect the seed function used in j2] 
which was of compact support. 


D. Solving the constraints 


where s l is the unit normal to the surface. Our approach 
to the apparent horizon search is based on that of [40] as 
also presented in mm- First given the spatial metric 
and extrinsic curvature / y$j, Kjj in Cartesian coordinates, 
we transform to work in spherical polar coordinates de¬ 
fined by 

r 2 = x 2 + y 2 + (z — zo ) 2 , 6 = arccos ^^ . 

(119) 


with 6 £ [0, 7 r] and where we take the z-axis to be the 
symmetry axis. In axisymmetry without twist, the spa¬ 
tial metric and extrinsic curvature then take the form 


Sij — 


( S rr 

r sin 0S r T 


V 0 


r sin 9S r T 

r 2 S 8T 

0 


: ) 

r 2 sin 2 dS^T ) 


( 120 ) 


Compactified coordinates: To solve for moment of 
time symmetry initial data, we write the spatial metric in 
spherical polar coordinates (r, 0, </>), and compactify the 
radial coordinate, leaving us with coordinates {A,Q, <j>). 
The compactification is defined either by, 


mA 

r= 2 (1 - A) ’ 

as suggested in [38], and used in [7] 
solver employed presently, or 


(116) 


in the same elliptic 


in the <f> = 0 plane. Local flatness on axis implies that 
the components S rr , S r T , S/t and S^t are even functions 
of 9 around the symmetry axis, with Sgx — S^t ~ 9 2 
around 9 = 0 , and similar dependence around 9 = n. 
Working in the p-z plane we may parametrize an appar¬ 
ent horizon by the level set s = 0 of 

s = r — F{9 ), (121) 

in terms of which the apparent horizon condition m 
can be rewritten as a first order ODE system, 


lA 


2(1-A 2 ) ’ 


(117) 


similar to that employed for example in |39j . The param¬ 
eter m partially controls the rate of compactification, but 
in either case spatial infinity corresponds to A = 1. 

Numerical solution: To discretize we employ a 
Chebyschev discretization in the radial A direction, and a 
Fourier grid in the angular directions. Since the Hamil¬ 
tonian constraint in this context is linear, solving the 
constraints amounts to a matrix inversion. With our 


particular method we find that the choice (116) leads to 


slightly worse constraint violations at a fixed resolution. 
One possible cause of this is that the coordinates (116) 


are irregular at the origin. Perhaps it is possible to use 
the alternative compactification in the construction of 
trumpet or puncture blackhole initial data, but we leave 
this for future consideration. 


E. Axisymmetric apparent horizons 

Formulation of the AH conditions: An apparent hori¬ 
zon is a closed two surface in the spatial slice, with unit 
outward pointing normal s l , with expansion, 

H = D i s i - K + sV A'y = 0, (118) 


F' = G, 

G' = (sin 2 0 7r 2 T - lrrl g T )F 2 L 2 qV(T k l0 D k s + LI< ZJ ). 

( 122 ) 


for F(9) and G(0), where the unit spatial vector s l and 
magnitude L are given by, 

s i = LD jS , L~ 2 = (D iS )(Djs) , (123) 


and qij = jij — SiSj is the induced metric in the level set. 
These expressions are evaluated in spherical polar coor¬ 
dinates. As noted elsewhere this parametrization is not 
completely general, only being sufficient if the apparent 
horizon is a ray-body containing the point zq. Regularity 
of an apparent horizon means that G(0) = G(tt) = 0 
Search strategy: Given the metric and extrinsic cur¬ 
vature we go about searching for an apparent horizon in 
the follow ing w ay. First we choose zo; r o and integrate 
the ODE (122) from 9 = 0 to 9 = 7r/2, with initial con¬ 
ditions -F(O) = ro and G(0) = 0. We simultaneously 
integrate backwards from 9 = n to 9 = n/2 taking as 
initial conditions F(ir) = ro and G(7r) = 0. If we have 
an apparent horizon the forwards (A + ,G + ) and back¬ 
wards ( F~ , G~ ) solutions will satisfy, 


A F = F + (n/2) - F~(tt/2) = 0, 
AG = G+(tt/2) - G~(tt/2) = 0. 


(124) 
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This gives a non-linear root finding task on the func¬ 
tion S : R 2 —► R 2 defined by 


S(z 0 ,r 0 ) = (AF,AG). 


(125) 


One complication is that the ODE system (122) needs 
to be regularized on the axis to impose our initial condi¬ 
tions. This is straightforwardly done by using the regu¬ 
larity conditions above, resulting in, 

p+ ' drier _ Ker \ p2 


G' = 


10T 


IrT 


2 2 ”"iv 


4'7rr ^yj'lrr) 


(126) 


at 9 = 0 and similarly at 9 = 7r. To arrive at this 
expression we have explicitly used the regularity condi¬ 
tion Sqt—S^t ~ 9 2 . In our numerical implementation we 
transform from Cartesian components, so this condition 
is automatically satisfied and we can instead use the con¬ 
dition in a slightly more complicated form involving 7 
and K^t- To the best of our knowledge this regulariza¬ 
tion of the coefficients has not been used before. The 
second step of our search is to iterate on (zo,ro) until 
we find a solution, or until the method fails. As an al¬ 
ternative strategy, it is normally proposed to integrate 
the ODE from 9 = 0 to 9 = tt then perform a bisection 
search on G(tt). We were unable to obtain satisfactory 
results this way because every surface except the appar¬ 
ent horizon itself diverges near 9 = tt, making the bisec¬ 
tion hopeless. Reasonable first guesses for zq would seem 
to be the position of the maximum of the Kretschmann 
scalar, or, if an apparent horizon was already found in a 
previous time-slice, the coordinate center of the previous 
horizon. 

Horizon mass: In twist-free axisymmetry the appar¬ 
ent horizon mass Mh is related to the area of the appar¬ 
ent horizon Ah as, 


M% = 


Ah 

167t 


(127) 


We can compute the area of the apparent horizon as a 
simple integral, 


Ah = 27r 


f L 1 y /ir 2 siw9 A9. 

Jo 


(128) 


where we have used the fact that apparent horizon is a 
surface of revolution. Here 7 is the determinant of the 
spatial metric in Cartesian coordinates. 

Simplifying assumptions: We are often interested in 
finding apparent horizons centered at the origin in space- 
times that are additionally reflection symmetric about 
the z = 0 plane. In this case we can trade our root¬ 
finding search above for a bisection search by simply fix¬ 
ing ^o = 0 and integrating the ODE (122) from 9 = 0 
to 9 = tt/2. Here we start the integration from dif¬ 
ferent initial radii r 0 until we find points about which 
which G(n/2) changes sign. We then bisect in ro to find 
the apparent horizon, where G(n/2) = 0. We typically 
choose the criterion G{ 7t/2) < 1(R 8 to end the search. As 
in the more general case, if we find many such surfaces 
we take the outermost as the apparent horizon. 


Numerical implementation: In practice we search for 
an apparent horizon as follows. During a bamps evolution 
we output the necessary components of the spatial met¬ 
ric and extrinsic curvature in the y = 0 plane at different 
coordinate times. For the integration of the ODE, we use 
various ODE integrators in the GSL [44] . To determine 
the apparent horizon accurately as fast as possible we use 
the explicit embedded Runge-Kutta Prince-Dormand (8, 
9) method, a high-order adaptive step integrator. When 
convergence testing we use a simple fourth order Runge- 
Kutta integrator. To evaluate the metric and extrinsic 
curvature at each point (r = F, 9) along the level set we 
use barycentric Lagrange interpolation inside each bamps 
sub-grid. For the root-finding we again use the GSL, now 
choosing one of the ‘hybrid’ algorithms that do not need 
the Jacobian of the system of equations being solved. 
In Fig. [3] we present the apparent horizon found using 
our method for a centered p = 0, amplitude A = 12, 
Brill wave initial data set, comparing it with that which 
we find using a stand-alone apparent horizon finder im¬ 
plemented in the MATLAB initial data code employed 
in gg. 


F. The analytic Cartoon method 

Here we discuss the implementation of the so-called 
Cartoon method g] for axisymmetry in a pseudospectral 
method for the Einstein equations. We assume that we 
are given the 3d system in a Cartesian coordinate system 
x l in which all variables are smooth, T £ C°°. The basic 
idea of the Cartoon method is to apply wherever possi¬ 
ble the same coordinates and discretization that lead to 
stable evolutions in 3d. Hence we compute the axisym- 
metrically reduced system in Cartesian coordinates and 
with Cartesian tensor components, without adapting co¬ 
ordinates and thereby avoiding the coordinate singularity 
at the axis. 

Concretely, the computational domain is chosen to be 
the cc-z-plane defined by y = 0. Partial derivatives d x and 
d z are computed as for the 3d system. What is missing 
are the points and the numerical data in the y-direction 
for the computation of d y . However, we can obtain the 
y-derivative by invoking axisymmetry, since the fields in 
the y = 0 , x-z-plane determine the fields for y 7 ^ 0 by the 
rotation symmetry. Similarly, it suffices to consider only 
the half-plane x > 0 and y = 0 while still using the same 
stencils for d x and d z as in 3d. 

The Cartoon method was first introduced for a Carte¬ 
sian BSSNOK pHrl - lTSl code using finite differencing [5]. 
The d y derivative was computed by adding ghost points 
in the y direction, so that identical 3d stencils could be 
used for 3d and axisymmetric 2d calculations. For a spec¬ 
tral collocation method, we could do the same and popu¬ 
late a 3d spectral element by rotation. There would still 
be significant gains in efficiency since only a 2 d subset of 
a 3d spectral grid consisting of many patches needs to be 
populated. However, it is also possible to derive analyt- 
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FIG. 3: In the left hand panel the apparent horizon for a centered A = 12 Brill wave, as found by our apparent horizon finder 
and a bespoke Brill-wave apparent horizon finder, are plotted. This data has been used as a standard test case elsewhere in 
the literature [401 143] . We compute the ADM mass as Madm = 4.67, which compares perfectly with Madm = 4.67 ± 0.01 
given in [3D]. The horizon mass is Mh = 4.66, again in agreement with the literature. In the right panel we show pointwise 
self-convergence labelled hy N = 25,100, 400 and 800, with N + 1 the lowest number of points in the series, and where we 
evolved with 2N + 1 and 41V + 1 to make the plot. Note that very few points are needed to show clean convergence because the 
surface varies slowly in 9. This also means that one can not reliably convergence test at high resolutions because the difference 
between the computed surfaces are essentially at the level of round-off. 


ical formulas for d y in terms of quantities in the y = 0 
plane only, so this is clearly the preferred way to proceed. 
To our knowledge this was first implemented in [491 . in 
that case for finite differences and the second order GHG 
system. For an arbitrary smooth tensor T, axisymmetry 
is given by the vanishing of its Lie derivative along the 
rotational vector, C^T = 0 . 

Off-axis, x 0. Let us consider various tensor types 
of interest, suppressing their t and z dependence. For a 
scalar, 


d y u(x, 0) = 0. 


(129) 


The second derivative does not vanish in general. For 
vectors and covectors (i / 0), 

d y v x (x, 0) = — v y (x, 0), dyV v (x, 0) = —v x (x, 0), 

d y w x (x, 0) = — w y (x, 0) , d y Wy(x, 0) = ~w x (x, 0) . 

(130) 

the derivative is equal to the components of the vector di¬ 
vided by radius, with x and y components interchanged. 
For a symmetric (0,2) tensor (say, the four-metric g a b ), 
at y = 0, x 0, 


dyQtt - 

— 0; 9 y g t , 

* = o, 

d y g zz = 

dyQtx 

= ~x 9tv ’ 

d y gt y = 

1 

— —gtx ^ 

X 

dyQxz 

1 

9yz j 
X 

d y g y z 

l 

— ~9xz 1 
X 

dyQxx 

2 

9xy> 

X 

®y9yy 

2 

9xy i 
X 

dydxy 

= — {9xx 

X 

9yy )• 



(131) 


Some components behave like scalars, some like covec¬ 
tors, and some show the two terms occurring in the Lie 
derivative, which may result in a factor two due to sym¬ 
metry. 

On-axis, x = 0. Axisymmetry by itself does not im¬ 
ply differentiability on the axis. Consider, for exam¬ 
ple, u(x, y) = p. We combine axisymmetry with the 
condition that in Cartesian coordinates T £ C°° in two 
ways. First, consider parity under (x,y) —> (—x,—y), 
which corresponds to a rotation by n around the 2 -axis. 
Because of axisymmetry, we have T(p, 0) = ±T(—p, 0) 
and d y T(p,0 ) = =pd v T(—p,Q). Since d y T is continu¬ 
ous, the limit p —> 0 exists. Hence for tensors that 
are even under this type of parity, the derivative van¬ 
ishes, d y T even ( 0,0) = 0. For tensors that are odd, the 
tensor vanishes, T 0 dd{ 0, 0) = 0, and d y T a dd{ 0, 0) is a reg¬ 
ular, finite value. We therefore impose that d y vanishes 
on the axis for even quantities and ask how we can com¬ 
pute the value for the odd quantities. 

From vanishing of the Lie derivative, we obtain rela¬ 
tions for the tensor components themselves, not for their 
derivative. For a scalar, there is no extra condition. Ex¬ 
amples for relations obtained from (130) (131) are, 


u*(0,0) = 0, Wi(0, 0) = 0 

5 ** ( 0 , 0 ) = g ty ( 0 , 0 ) = g xz ( 0 , 0 ) = ^( 0 , 0 ) = 0 , 
g xy (0, 0 ) = 0 , g xx { 0 , 0 ) = g yy ( 0 , 0 ). ( 132 ) 

Although we obtain some of the same information that 
we already discussed for (x, y) -A (—x, —y ) parity, for 
even parity quantities with two or more indices there are 
additional relations. For the metric components these are 
related to covariance under rotation by 7r/2, or (x,y) — I 

To find the deriva tive d,, at ( 0 , 0 ), we invoke l’Hopital’s 
rule. Basically, in ( 130 ) ( 131 ) the - factors become a 
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partial derivative in x because the other terms vanish. 
For example, 

d y v x ( 0 , 0 ) = -d x v y (0, 0 ), d y v v ( 0 , 0 ) = d x v x (0,0). 

(133) 

Notice that starting with two-index components this is 
not entirely trivial since there is more than just one term 
on the right-hand side. 

Axisymmetry for partial derivatives of tensors. There 
also are objects like &i a b = dig a b, which are not ten¬ 
sors, but partial derivatives of tensors. The Lie derivative 
£- 4 >dig a b is in general not defined for non-tensors, and a 
priori it is not clear whether C.^dig ab = 0 implies axisym¬ 
metry. However, we can obtain the required formulas by 
computing 

SiG^g a b — £~ , (f)digab A gcb^a^if 1 T Qac&bdif) > (134) 

where C^ is introduced to collect those terms that corre¬ 
spond to the Lie derivative of a tensor, and the remaining 
terms are the deviation from the tensor formula. Note 
how the last term in d z (4> c d c g ab ) = (t> c d c dig ab + d c g ab di4i c 
provides precisely the term that would otherwise be miss¬ 
ing in the sum over index locations in C,^dig ab . 

The key observation is that in the case of a rigid 
rotation in adapted coordinates generated by cj> a = 
( 0 , — y, x , 0 ) T , the second derivatives of (j) a vanish, 


d a d b (t> c = 0. (135) 


Therefore, in this special case we obtain the correct result 
using the tensor formula, 

(di^'^gdb — ^(fi^igab^ (136) 


as was also noted in [5]. This generalizes immediately to 
partial derivatives of arbitrary tensors, and also includes 
the case of the Christoffel symbol requ ired for the BSS- 
NOK or Z4c system, compare 0. Eqn. (1351 furthermore 
simplifies the computation of second derivatives. 


VI. CODE VALIDATION 

In this section we present a set of numerical exper¬ 
iments performed to try and obtain an optimal setup 
for the first order generalized harmonic system for our 
gravitational wave collapse evolutions that follow in later 
work. 


In the following set of experiments we always take A = 
0.01 and cr = 10, and fix the grid setup. We take the stan¬ 
dard formulation used in the SpEC code, namely 70 = 
—71 = 72 = 1, and 74 = 75 = 0. We impose outer 
boundary conditions at a coordinate radius of r = 16 
and evolve in 3d with octant symmetry imposed. 

Harmonic gauge: Starting with the pure harmonic 
gauge H a = 0, we find that the outgoing gauge wave 
is harmlessly absorbed using either the gauge boundary 
condition (31) or (32). At the particular resolution and 
grid-setup that we chose for this test the harmonic con¬ 
straint violation at the end of the evolution, t = 100 , is 
around 10 -14 and shows no sign of increasing with ei¬ 
ther choice of gauge boundary condition. The difference 
between the results with the two gauge boundary condi¬ 
tions is rather small, the maximum difference in the shift 
being around 10~ ‘ at the end of the run. But here the 
initial pulse is very weak, and this is of no concern. In the 
left panel of Fig. [4]we plot | a — 1| in the outer boundary, 
to demonstrate how the coordinates settle down. 

Generalized harmonic gauge: Switching now to use 
the generalized harmonic gauge condition (101 with = 
0.4, p = 1 and 77 s = 6. Using then the gauge boundary 
condition (31) we find that the shift starts to grow at the 
boundary, and the numerics fail at t ~ 42. This behavior 
is perhaps not surprising given the large damping coeffi- 
cents and the understanding obtained for the constraint 
preserving subsystem with damping in section IIB The 
gauge source functions have the same effect on the gauge 
as the damping terms on the constraints, namely they 
cause reflections from the boundary. We expect that it 
will be suppressed as the outer boundary is placed fur¬ 
ther out so that the gauge sources are smaller where the 
boundary condition is applied. Using instead the gauge 
boundary conditions (32) this growth is completely ab¬ 
sent, which is why we do not implement conditions de¬ 
rived explicitly to reduce gauge reflections in the present 
work. This behavior is demonstrated in the right panel 
of Fig. [4] where one sees the magnitude of the shift vec¬ 
tor in the outer boundary in each case. With the gauge 
boundaries (32), at the end of the run the harmonic con¬ 
straint violation C x is around 10~ 14 and appears not to 
be growing. Looking at the shift however, it does seem 
that some further improvement may be possible in the fu¬ 
ture, as its peak lies at the outer boundary, with a value 
around 10 -11 . 


B. Constraint experiments 


A. Gauge boundary 

Gauge wave initial data: We evolve the Minkowski 
line-element with a perturbation initially placed in the 
lapse, so that, 


(137) 


Simplified subsystem: We now repeat some of the ex¬ 
periments of the previous section with the choice 74 = 
75 = 1 / 2 , and with different choices of 70 , using al¬ 
ways the gauge boundary condition (32). With the 
pure harmonic gauge H a = 0, we find that the con¬ 
straint violation at t = 100 is again around 10“ 14 if we 
take 70 = 1 , and slightly larger, but still less than 10“ 13 
if we choose 70 = 0 . 02 , the value suggested by the exper- 


a(t = 0) = 1 + A e -[x 2 +y 2 +z 2 ]l° _ 
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FIG. 4: In the left panel we plot the |a — 1| in the outer boundary as a function of time, obtained in the evolution of a gauge 
pulse on flat space, initially centered at the origin. The coordinates eventually seem to settle on, or very close to Minkowski 
slices. On the right we plot the magnitude of the shift in the out er b ou ndar y using the harmonic damped wave gauge to evolve 
the same gauge pulse with either the gauge boundary condition (311 or (32 I. In the former case the shift rapidly grows, causing 
the code to crash. 


iments in [50] for a related formulation. Moving to the 
generalized harmonic choice (101 once more, we find that 
again that the violation at the end of the experiment is of 
the same order as when using the pure harmonic gauge. 
The result is plotted in Fig. [5] These results may not be 
representative when evolving different initial data, but 
we cautiously take 74 = 75 = 1/2 and 70 = 1 as our de¬ 
fault setting, periodically testing different choices, most 
often playing with 70 in such experiments. 

Constraint preserving conditions: We performed the 
same experiments, with the generalized harmonic gauge 
and the new default formulation parameters, changing 
to the alternative constraint boundaries (39) or (451 and 
found first that the violation throughout is very similar 
to the initial choice (28 1 . Although initially the violation 
with the reflection reducing condition is slightly smaller 
than with the ‘geometric’ condition, later on there is 
practically nothing to choose between them. Consider¬ 
ing that the violations are in the round-off regime 10 -14 
it is hard to judge from this experiment which of the 
conditions behaves most favorably. 


C. Lapse power in constraint damping 


Initial data: We now evolve centered A = 2.5 Brill 
wave initial data, which is subcritical, with an ADM 
mass of Madm = 0.19. We evolve on the same grids 
used in the previous section, but with a slightly higher 
resolution (19 3 rather than 15 3 points per cube). We 
evolve using 70 = 0 . 2 a* with l = 0 , the standard choice 
elsewhere, or l = — 1 , a modification which we hope will 
reduce constraint growth in the strongest field region. As 
above we use the generalized harmonic gauge (10 1 . We 
use only the gauge boundary condition (321. 


Basic dynamics: The Kretschmann scalar initially 
has a peak at the origin, evaluated around 2300 on the 
bamps grid, slightly less than in the previous study I7J. 



FIG. 5: We show the C x component of the harmonic con¬ 
straint along the x at time t = 100 for two different sets of 
constraint damping parameters with formulation parameters 
74 = 75 = 1 / 2 . in the evolution of a gauge pulse on flat-space 
as in Fig. [4] with the generalized harmonic gauge. On this 
basis we take these formulation parameters with 70 = 1 as 
our standard choice. 


This peak oscillates at the origin, peaking after an initial 
bounce with value around 500. The feature then rapidly 
propagates away and by a coordinate time t = 10 , the 
peak value on the grid is less around 10 -2 . The lapse 
initially decreases at the origin, this feature then propa¬ 
gating out to the outer boundary, behind which the lapse 
drifts back towards its initial value, unity. 

Constraint violation: Examining the C x constraint 
for the A = 2.5 data along the ir-axis, we see that only 
very small differences in the constraint violation between 
the l = 0 and l = — 1 evolutions. The small differences are 
not surprising because the lowest value the lapse function 
takes is around 0.78 having started from 1. The peaks of 
the C x constraint in the l = — 1 evolution are about 2-5% 
smaller than in the l = 0 run. Increasing the amplitude 
of the initial data to A = 4, one might expect the im¬ 
provement to be more significant as the lowest value of 
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FIG. 6: Comparison of the results of a Brill wave A = 1 evolu¬ 
tion with BAM and bamps. We show snapshots of the metric 
component along the x axis at t — 1.625. In the upper 
panel we show the pure harmonic gauge, and underneath the 
damped wave gauge with t/l = 0 and r/s = 1.0. The results 
of the codes are in good agreement in either case. 


lapse decreases to 0.37, but the difference still amounts 
to between 2-5% at the peaks of the violation. 


D. BAM vs. bamps comparison 


Another validation strategy for bamps is to compare 
the numerical results with those of an independent code. 
For this we used BAM ED, evolving identical initial data 
with the same gauge conditions. This comparison we 
performed by evolving a centered Zq = 0 Brill wave 
with A = 1. We chose this weak amplitude because 
evolving the Brill data accurately with BAM rapidly 
becomes expensive as A increases in magnitude. We 
used pure harmonic slicing tjl = 0 with either harmonic 
shift r]g = 0 or the damped harmonic shift r)s = 1. In 
the BAM code we evolve with the BSSNOK formulation, 
for completeness, this gauge condition is given by 


dtp = a 2 \ f* + \P j dj ln X - 7% In a 

-nsp + pdjp. 


(138) 


in terms of the conformally decomposed BSSNOK vari¬ 
ables. For this test we did not employ the spherical shells 
or constraint preserving boundary conditions of 152] . 
Since the outer boundaries were placed at x = y = z = 
12, the solutions to the continuum PDEs being solved are 
not identical. Therefore we should not hope for perfect 
agreement for long. In Fig. [6] we plot the spatial metric 
component j xx at t = 1.625, when the agreement is still 
very good for either choice of the shift, being practically 
indistinguishable by eye. In practice the main source of 
disagreement at the resolution of these runs comes from 
mesh-refinement boundaries in the BAM grid setup. 


E. Octant and Cartoon 


Initial data and grids: To test our implementation 
of symmetry reduced expressions, either octant, Car¬ 
toon, or their combination we evolve weak A = 1 cen¬ 
tered pure plus polarization initial data as described in 
section [VBj using once again the generalized harmonic 
gauge (10) and the gauge boundary condition 32 We 


started with a base cubed sphere 3d grid with N = 15 
points per direction, and the number of subpatches de¬ 
rived from Afcu = 5, J\f cs = 4 and J\f cs = 3. The outer 
boundary was placed at r = 12 in the units of the code. 
For ease of comparison, the breakdown of the grids was: 



^y"total 

^y"total 

^"total 

^ytotal 

^/ytotal 

3d 

125 

600 

450 

1175 

4 x 10 6 

octant 

27 (12,6,1) 

48 (48,12) 

81(36,9) 

216 

5 x 10 5 

Cartoon 

25 

80 

60 

165 

4 x 10 4 

cart. oct. 

9(4,1) 

24 (8) 

18(6) 

51 

10 4 


where the numbers in parentheses denote the number of 
those grids that were cut in half (at the axis) once, twice, 
or three times respectively, for the 3d grids, and once 
or twice for the Cartoon grids. Note that our current 
non-octant Cartoon implementation is not optimal be¬ 
cause we evolve the whole x-z plane, wasting effectively 
a factor of two. Currently we use the code most often 
in Cartoon octant mode, so fixing this does not have a 
high priority. Looking at the table the main observation 
is that the expected reduction factor of eight (four) in 
the total number of grid points is present between the 
3d (Cartoon) and octant grids, but that this number is 
not so closely reflected in the grid breakdown, where we 
get only a factor six (three) in the total number of grids. 
This is obviously because there are many grids with fewer 
points. Since our parallelization does not take this fact 
into account, it is possible that one MPI process is given 
all non-cut grids, and so we can expect that the speedup 
rate is determined to a large extent by ratio in the num¬ 
ber of grids. As we make the domain larger the relative 
number of cut grids decreases, so we might expect that 
asymptotically the full speedup factors of eight or four 
can are attained. 

Basic dynamics: Although irrelevant for the octant 
Cartoon comparison, since these data have not been 
used before, we give a brief description of their evolu¬ 
tion. Initially the peak of the Kretschmann scalar occurs 
at p = ±0.65 with a value 1.18. This profile then oscil¬ 
lates about three times at the origin, attaining a peak 
value of 7.25 before rapidly dispersing. Looking at the 
lapse we see the familiar behavior that at the origin it 
oscillates slightly before presenting a longer decrease, al¬ 
though at the minimum is only 0.995, having started 
from a(t = 0) = 1 everywhere. Afterwards this pulse 
propagates out, roughly following the disturbance in the 
Kretschmann. Looking at the shift component /3 X along 
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the a;-axis we find that early on there is a growth which 
peaks at x = 1.06, with value 0.0027. The development 
of the shift looks more like a slowly oscillating standing 
wave than a localized propagating feature. 


3d, octant, Cartoon and octant-Cartoon comparison: 
Taking first the 3d and octant evolutions, we see near 
perfect agreement throughout the evolution. There are 
small differences however, starting from the beginning 
of the simulation at the level of round-off; differences 
of 10 -15 in metric components, which slowly drift as the 
evolution goes on. This behavior is expected because 
the derivative approximation differ at this level. Similar 
differences were found between the other setups. These 
differences are never larger than the constraint violation, 
in for example C x , and we have looked at convergence 
(see section VIF for more discussion) with each setup, 
although not for this data, and find no indication of a 
problem. For the speed comparison we ran the code with 
each setup on 24 cores (with hyperthreading) of our local 
cluster Corel2 with Intel Xeon X5650 processors. The 
octant run was a little more than 6 times faster than 
the 3d run, as expected given the foregoing discussion. 
The octant Cartoon run was about 2.4 times faster than 
the pure Cartoon test, which is a little disappointing. 
Going from J\f ss = 3 to J\f ss = 6 radial subdivisions in the 
outer shells, this value increases to 2.9, demonstrating 
the expected dependence. Comparing the full 3d and 
octant Cartoon runs, there was a gratifying speed up of 
nearly a factor 400. 


F. Convergence 


The bamps numerical method gives us two options for 
increasing resolution. The first is to add grid-points in 
every domain, the second is to subdivide grids further, 
keeping the number of points inside each subpatch fixed. 
Given fixed finite computational resources it is not obvi¬ 
ous what is the optimum strategy to achieve the smallest 
possible error, because although we would expect that 
adding points brings spectral convergence, it also comes 
with a N~ 2 dependence in the allowed time-step, whereas 
on the other hand, as we will see, adding more sub¬ 
patches allows the code to scale up to a large number 
of processors. Probably the optimal strategy relies on 
a balance between each. To examine the effect of each 
strategy in the simplest possible way, we evolved gauge 
wave initial data on the Minkowski spacetime, which 
was setup by choosing a = 1 + A exp[— (r/cr) 2 ], (3 l = 0, 
with r = \/x 1 + y 2 + z 2 as usual, and otherwise the flat 
spatial Cartesian metric and vanishing extrinsic curva¬ 
ture. The results are plotted in the four panels of Fig. [7] 
and confirm our expectations. 


G. Filtering 


To demonstrate the necessity of the filter (561 we 


evolved a centered A = 1 Brill wave. The results are 
plotted in Fig. [8] In the left panel we see that without 
filtering the constraint violation starts to grow exponen¬ 
tially in time, whereas with filter the growth is completely 
absent and the norm of constraints remains steady at a 
very low value. In the right panel we plot the magni¬ 
tude of the fourth highest spectral coefficient of g xx in 
the transition shell as a function of time. This coefficient 
is the first that is directly unaffected by the filter. We 
see that the growth in the constraints seems to be associ¬ 
ated with an explosion in the higher spectral coefficients. 
Interestingly we tried the same experiment with gauge 
wave initial data, and did not see the effect, at least in 
the same time-frame. We expect that the same behavior 
would manifest if we were to evolve long enough. The 
obvious conclusion we draw from this is that it is impor¬ 
tant to test these methods with several data types to get 
a reliable picture of their properties. 


H. Performance 

Strong-scaling: The current bamps parallelization 
strategy is to obtain perfect scaling using many sub¬ 
patches, and splitting these subpatches across many pro¬ 
cessors. The key is that, contrast to buffer zones required 
in the decomposition of a finite differencing grid, only 2d 
surfaces of points need be passed by network communi¬ 
cation, making the relative time spent there negligible. 
In a finite differencing approach the relative size of the 
buffer zones decreases with resolution, but in practice 
can still be significant in production runs. In Fig. [9] we 
present strong scaling plots performed on the SuperMUC 
cluster located in LRZ Garching, with Intel Xeon E5- 
2680 8C processors. We ran the code in 3d. We took a 
grid with 4459 total subpatches, and increased the num¬ 
ber of cores used until we were computing one patch per 
core. We find perfect scaling. On the other hand bamps 
is currently not parallelized whatsoever at the subpatch 
level, which means that the maximum number of points 
per subpatch is in principle determined by the amount of 
memory available to one core. At least when running the 
code in Cartoon mode however we do not find, in prac¬ 
tical terms, that this is problematic. Instead the N ~ 2 
restriction in the time step makes increasing the number 
of points infeasible long before we are close to filling the 
available memory. In 3d this may no longer be the case. 
We leave such considerations for future work. 


VII. SINGLE BLACKHOLES 

The main thrust of our development has been towards 
treating collapsing axisymmetric gravitational waves ac¬ 
curately. For super-critical data the cubed-ball grid is 
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FIG. 7: Evolution of a gauge wave with A = 0.01 and a = 1.0. In the upper panels we used a spatial resolution of IV = 21 on 
a grid with M = 1 subpatches. The upper left panel gives a snapshot of gtt along the x axes at t=3.55. The upper right shows 
the Chebyshev expansion coefficients at the same time with the same color coding. The lower panels show convergence of the 
constraints for the same initial data; on the left we increase the number of points N in each grid, on the right we increase the 
number of subpatches A f. 


unsuitable after the formation of an apparent horizon. 
Therefore the strategy for long-term evolution is to take 
the data and interpolate onto a cubed-shell grid, with the 
excision surface suitably positioned, changing the lapse 
and shift to be sure that the excision surface is a true 
outflow boundary. A necessary requirement is to treat a 
single blacklrole, which is what we discuss here. 


curvature take the form, 


g ab dx a dx b 


2 M\ l2 AM 

- df 2 +- dt dr 

r J r 



dr 2 + r 2 clfl 2 , 


(139) 


with dfl 2 the flat metric on the two-sphere, and 


A. Initial data 


Kjjdx'dx- 1 


2 M 




(140) 


Kerr-Schild coordinates: We evolve the Schwarzschild 
solution in Kerr-Schild coordinates as was done with an 
earlier version [M] of the present code. Although the cur¬ 
rent numerical method is not particularly close to that 
used previously, some components of the older code were 
inherited. Importantly evolving this data allows a sim¬ 
ple comparison with the previous method and results. 
In spherical polar coordinates the metric and extrinsic 


respectively. Inside the code the line-element is written in 
Cartesian coordinates in the standard way. More discus¬ 
sion of Kerr-Schild coordinates can be found in [Ml 155] . 

Harmonic Killing coordinates: We additionally 
evolve starting from the harmonic Killing slicing de¬ 
scribed in [56], which serves as a convenient starting 
point when transitioning from one generalized harmonic 
gauge to another. For this initial data, in spherical polar 
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FIG. 8: Influence of the filter at example of a A = 1 Brill wave evolution. On the left we show the time evolution of the 
constraint monitor C mo n- In the simulation using a filter the constraint violation settle down to 10 -1 °. Without using a filter 
the constraint violation grows and lead to a failure of the simulation at t ~ 150. On the right we show the evolution of the 
fourth highest Chebyshev expansion coefficient. It is the highest mode which is not affected by the filter. Without the filter 
the high frequencies grow over time and cause the simulation to fail. The filter sets the highest frequency to zero which avoids 
the growth of the high frequency modes. 



FIG. 9: Strong scaling of bamps with Af = 5 on the Super- 
MUC cluster. Here a grid with Af = 5 sub patches was used. 
In total this grid consists of 4459 patches. 


Cartesians according to, 

x = (r — M) sin 9 cos 4 >, y = (r — M) sin 9 sin (j >, 
z = (r — M) cos 9 , (143) 

The resulting metric has a coordinate singularity at r = 
M, with r implicitly defined in the obvious way from 
the new coordinates. The coordinate singularity is not 
a principle problem as we could just put the excision 
surface outside this radius. But bamps relies on standard 
Cartesian coordinates in several places. So in the code 
we could transform in the standard way but then choose 
the gauge source function, 

H a = 2{JdJf ab \. (144) 


coordinates, the metric and extrinsic curvature are 

g ab dx a dx b = - (l - —) d t 2 + ^ dt dr 
V r J r z 

, 4 M 2 \ f 2 M\ <2 2,^2 

+ ( H- 2 — I 1 1 +- I dr 2 + ? ,2 dfl 2 


(141) 


and 


K = — 


Kgg = 


4 M 2 4M 3 + dM 2 r + 3 Mr 2 + 2 r 3 


1 + ^/1 + ^ 


1 + 2M U + 4M1 


(142) 


with the remaining components vanishing. For this data 
spatially harmonic coordinates are obtained by building 


with J a a / the Jacobian between the standard a-index 
Cartesians and harmonic Cartesian a! index coordi¬ 


nates (143), the compound object ( JdJ) is defined by, 


(JdJ) a bc = (J- 1 )? dj a a ,. 


(145) 


with J“, = and where indices are manipulated 

in the obvious way with g a g to obtain (144). Instead we 


just choose the gauge source function to be fixed at its 
initial value, as will momentarily be discussed. In this 
section we use the code exclusively in Cartoon mode, on 
a cubed sphere grid. We start with the excision surface 
at r = 1.8 M, and the outer boundary at r = 31.8 M. 
In our base setup we take Af = 3 radial subpatches each 
with N = 25 points per direction. The runs were per¬ 
formed on a desktop machine with an eight-core intel i7 
CPU, which was able to compute at about 250M/hour, 
the base run requiring about 14 MB of RAM. 
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B. Freezing gauge source functions 

Killing gauge sources: Given initial data which ad¬ 
mit a time-like Killing vector, we can ensure that the 
evolution of the system is trivial, at the continuum level, 
neglecting the effect of outer boundary conditions, by 
choosing the Killing lapse and shift, and taking the gauge 
source functions H a so that d t a = <9 t /3 z = 0 initially. In 
particular we must choose, 

H a = -T a (t = 0), d t H a = 0. (146) 

Kerr-Schild evolutions with SpEC GHG: We began 
by evolving the Kerr-Schild initial data with the stan¬ 
dard formulation parameters of 0 , namely 74 = 75 = 0 
and 70 = 1 on our base grid as just described, using 
the gauge boundary conditions ©. Immediately we see 
that the innermost subpatch has the largest constraint vi¬ 
olation, peaked at around 1CP 6 in the C x component of 
the harmonic constraint. This is not surprising because 
the innermost subpatch contains the part of the solution 
with the largest derivatives. The evolution successfully 
continues until the final time t = 1000 M. But after the 
initial expansion to 10~ 6 , a slow expansion in C x is visi¬ 
ble, and this growth becomes more rapid as the simula¬ 
tion continues. By the end, the maximum value of C x is 
around 10^ 3 , with peaks appearing at the inner and outer 
boundary of roughly the same size. We then increased 
resolution from the base grid to N = 27, 29 and N = 31. 
The N = 27 point grid runs at about 178 M/hour, and 
the initial peak in the C x constraint violation is reduced 
by a factor of about two, with this ratio of improve¬ 
ment slowly declining until the end of the evolution. 
The N — 29 grid runs at 129 M/hour, with both the 
initial magnitude of the violation and the ‘slow expan¬ 
sion’ of the C x constraint quashed, the peak being a 
factor 2.8 smaller than in the base run at the end of 
the simulation. The highest resolution N = 31 point 
grid ran at 96 M/hour, with the final improvement in C x 
against the base run being a factor of 5.3. Since the 
largest constraint violation occurs in the excision sub¬ 
patch an obvious question is whether or not the excision 
and outer boundaries would interact badly if they were 
on the same grid. Although the issue is of little practical 
concern for production runs, for development it deserves 
a little attention, and therefore we evolved our base grid 
from before, but cutting the outer two subpatches so that 
the outer boundary lies at 11.8 M. This test is not com¬ 
pletely fair because the outer boundary conditions are 
expected to perforin better as they are applied further 
out. We find that the initial peak in the violation of 
the C x constraint is about five times greater than in the 
base run at t = 200 Af. At the end of the evolution again 
at t = 1000 M by coincidence the constraint violation in 
the restricted domain is smaller, but this is just because 
the slow oscillations in each simulation are out of phase. 

Kerr-Schild incoming wave evolutions with SpEC 
GHG: Next we evolved the same initial data and gauge, 
but this time with the same domain as in Fig. 3 of [8]. To 


do this we took Af = 2 radial subpatches, with the same 
base resolution as before, so that the outer boundary is 
placed at r = 21.8 Af. We similarly specify exactly the 
same given data for an incoming gravitational wave as in 
that study, taking in particular, 

d t h ab = f(t)(x a x b + y a y b - 2z a z b ), (147) 


with the vectors here the coordinate vectors defined in 
the obvious way. We take, 


f{t) = Aexp[—(f - t P Y/u ) 2 }, 


(148) 


10 


we 


with A = 10 , t p = 60 Af and u> = 10 Af. In Fig. 

show the results from these experiments, obtained with 
a sequence of different resolutions. We plot the Weyl 
scalar ’F 4 (24), averaged over the outer boundary, 


4tt (f?^ 4 ) 2 


|^ 4 | 2 dA. 


(149) 


The surface area of the outer boundary is 4ttR 2 . Fitting 
the highest resolution data between t = 100 and t = 
200 we find a ring-down frequency of 3?[wAf] ~ 0.372 
as expected [57j. In this evolution we found that the 
apparent horizon oscillates slightly as the gravitational 
wave is absorbed, increasing the horizon mass ( 127| by 
about 6 x 10 - 7 M , with Af the ADM mass of the analytic 
initial data. Note that the gauge boundary condition 
being employed here is not identical to that used in [ 8 ], 
so the agreement is remarkable. The effect of much larger 
pulses of gravitational radiation falling onto a blackhole 
using similar gauge conditions was studied in [581 . 

Discussion of and comparison with m The prior 
bamps study focussed on obtaining numerical stability 
in the evolution of a single Schwarzschild blackhole with 
the Kerr-Schild slicing. The numerical method used a 
Chebyscliev-Fourier-Fourier spatial discretization on a 
single shell with a spin weighted spherical harmonic fil¬ 
ter to prevent high frequency growth of the error. In 
that study the outer boundary condition employed sim¬ 
ply fixed the incoming characteristic variables (| 6 ]) to some 
given data, namely their initial values. This approach 
is possible only when the analytic solution is known, 
otherwise incoming constraint violations are generated. 
Placing the inner boundary at r = 1.8 Af and the outer 
boundary at r = 11.8 Af, very long evolutions, say until 
at least t = 200 000 Af, could be performed with little res¬ 
olution, in accordance with [ 8 ] . On the other hand, using 
this method, the naive boundary conditions rapidly dete¬ 
riorated as the outer boundary was pushed out, and, cru¬ 
cially resolution did not help but rather made the prob¬ 
lem worse. A possible explanation for the latter effect is 
that no filter was being applied in the radial (Chebyschev 
discretized) direction, which have already seen is a cru¬ 
cial ingredient for stability with the current method. The 
likely cause of the boundary problem is that, as explained 
in nn. boundary conditions that just freeze the incom¬ 
ing GHG characteristic variables are orders of magnitude 
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FIG. 10: The right panel shows the average over the 
Weyl scalar 'hi in the outer boundary in the evolution of 
Schwarzschild perturbed by a small gravitational wave in¬ 
jected through the boundary. In the left panel we see con¬ 
vergence of the constraints as resolution is increased. At 
lower resolutions a drift is present in the ring-down. There 
is good agreement with Fig. 3 of [8], and the ring-down fre¬ 
quency agrees well with the analytical computation [57]. At 
the end of the test there is some disagreement with [8], but 
since square-roots of very small quantities are being taken 
we expect this is caused by round-off error. It seems that 
on the cubed-sphere grid more resolution is needed to obtain 
clean results than with the spherical harmonic discretization 
used in [8]. This is perhaps not surprising, since the latter 
discretization is well-suited to the given data. 




more reflecting than the Sommerfeld like choice contained 
in (31). Evidence for this is obtained in the current code 
by changing from the gauge boundary condition (31) to 
use instead, 


^ )cd [dtu~ d ]= 0, (150) 

evolving once more the Kerr-Schild initial data on the 
base grid. Placing the outer boundary further out then 
results in greater reflections. However rather than trying 
to improve a condition only suitable for evolving known 
data, we immediately moved to the constraint preserv¬ 
ing, radiation controlling conditions, with which this is¬ 
sue is completely absent. The first attempted implemen¬ 
tation of a regular center in the bamps code was to use the 
Chebyschev-Fourier-Fourier discretization with a double 
covering in the radial direction, similar to that employed 
in [59]. The approach was not successful, as we always 
eventually found irregularities in the numerical solution 
at the origin. An exponential filter was applied to the 
Chebyscliev coefficients in the radial direction, but to 
little effect. Eventually we settled on the cubed sphere 
approach, in part because of the expectation that they 
will later be more convenient for mesh-refinement. Other 
possible solutions to the problems we faced would be to 
use one-sided Jacobi polynomials as in SpEC [BO] or to 
employ a filter that projects the solution in another basis 
onto these polynomials. 





FIG. 11: Comparison of the evolution of Schwarzschild 
with Killing-Kerr-Schild gauge sources with either the gauge 
boundary condition (311 or the alternative (32 I at the end of 
the simulation t = 1000 M. In the upper panel we plot the 
logarithm of the constraint violation C x - In the latter case 
the violation is greatly reduced. In the lower two panels we 
show the lapse and shift; the drift present when using ( |31| ) is 
practically absent with (32 I. 


Kerr-Schild evolutions with simplified constraint sub¬ 
system: Using our standard choice for the formulation 
parameters 74 = 75 = 1 / 2 , and taking 7 0 = 0 . 2 , return¬ 
ing to our base resolution from the tests with the SpEC 
version of GHG, we find that by t = 200 the C x constraint 
is about 5 times larger than that we obtained before, and 
by the end of the simulation the new run has accrued 
a C x constraint violation with a sharp peak at the outer 
boundary of order 10 _1 . This result seems to be in con¬ 
tradiction to those of section |VIB[ until we remember 
that there the gauge boundary condition (321 was em¬ 
ployed. Increasing the constraint damping to 70 = 1, the 
initial violation is comparable to the SpEC GHG evolu¬ 
tion previously described throughout the evolution, and 
the spike at the outer boundary is suppressed by roughly 
an order of magnitude. At the end of this run the maxi¬ 
mum of the C x constraint occurs at the excision bound¬ 
ary with a value around IIP 3 . This experiment thus 
highlights that the choice of the damping parameters and 
boundary conditions can be rather subtle. 
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Kerr-Schild evolutions with alternative boundary con¬ 
ditions: Next we returned to the base grid, and switched 
to the alternative gauge boundary conditions (32 1 , 
with 74 =75 = 1/2 and 70 = 1. We find that the 
aforementioned growth in the constraints is completely 
eradicated, and the drift in the lapse and shift is also 
suppressed. Evolving the same data with the same for¬ 
mulation and gauge boundary condition, but using the 
modified constraint preserving boundary condition (39) 
gives almost identical results. Using instead the reflection 
reducing conditions (45) we see a small improvement in 
the violation throughout the simulation. Repeating the 
experiment with the incoming gravitational wave injected 
through the boundary with the standard constraint pre¬ 
serving condition (28) and the gauge boundary condi¬ 
tions (32), the growth visible in Fig. 10 is also completely 
absent, even on the base resolution N = 25 grid. 

Harmonic Killing slice evolutions: We now returned 
to our base grid and resolution, taking the formulation 
parameters 74 = 75 = 1/2, and 70 = 1, evolving the 
Harmonic Killing slice with the gauge boundary condi¬ 
tion (31). The test successfully runs to t = 1000 M. 
Comparing with the equivalent evolution of Kerr-Schild 
data, we see that initially near the excision boundary 
the C x constraint violation is significantly greater in the 
Harmonic Killing test. By t = 200 M this difference has 
accrued to around two orders of magnitude. Later how¬ 
ever, as the violation in the Kerr-Schild Killing evolu¬ 
tion starts to grow, it overtakes that of the Harmonic 
Killing evolution. At t = 1000 M the peak of the con¬ 
straint violation in the Harmonic Killing run is about an 
order of magnitude smaller than in the earlier test. As 
remarked before, in the Kerr-Schild test the inner and 
outer boundaries have roughly the same magnitude in 
the C x constraint violation. Interestingly the twin peaks 
are not present in the Harmonic Killing data because 
the outer boundary is hugely improved. This finding is 
consistent with the gauge wave tests presented in sec¬ 
tion [VIA] although this test is somewhat easier for the 
gauge boundary conditions because of the complete lack 
of dynamics present in the gauge wave test. In the Har¬ 
monic Killing evolution we are evolving with pure har¬ 
monic slicing, and some non-zero spatial gauge source 
functions, which suggests perhaps that the growth at the 
outer boundary is predominantly caused by the use of a 
non-trivial gauge source function for the lapse function, 
as it interacts with the boundary. Indeed looking once 
more at the lapse function towards the end of the Kerr- 
Schild evolution we see that it is drifting from its initial 
value, but that this effect converges away with resolu¬ 
tion. In any case, the peak in the constraint violation 
at the outer boundary in the Killing Kerr-Schild data is 
suppressed as the outer boundary is placed further out. 


Harmonic Killing slice with gauge perturbation: A de¬ 
sirable property for a set of dynamical coordinates is 
that in the presence of a, perhaps approximate, time¬ 
like Killing vector they quickly asymptote to a time- 
independent state. For an arbitrary physical or gauge 


perturbation there is no hope that this will occur, and 
nor can any finite set of numerical experiments prove 
that that there is a basin of attraction to a stationary 
state. We can however look for some indication of this 
behavior. To do so we start by taking the initial data 
for the Killing harmonic coordinates, and then perturb 
the initial lapse function by Gaussian as in the previous 
gauge wave evolutions. In terms of the first order GHG 
variables this is a slightly fiddly procedure, as compared 
with the use of lapse, shift and spatial metric, so we give 
a quick summary: 

• Set spatial metric and extrinsic curvature from the 
exact solution. 

• Take the Killing lapse and shift. Use the condi¬ 
tions d t a = 0 and d t ff = 0 to set the gauge source 
functions H a . 

• Add the desired perturbation to the lapse (or shift) 
and then transform to the first order GHG vari¬ 
ables. 


We perturbed the lapse by a Gaussian, 
Aa = A exp [ — 2 (r — ro) 2 ] 


(151) 


with A = 0.3 M and ro = 4 M. A similar experiment 
was made in }61j . but starting from a Maximal slice 
of Schwarzschild to test the gauge driver system. We 
find that the perturbation in the lapse propagates away, 
rapidly leaving behind the solution with the harmonic 
Killing data with unperturbed spatial coordinates, or at 
least negligibly perturbed. The greatest danger to the 
evolution is probably that the excision boundary fails to 
be outflow, but at least with this perturbation that does 
not occur. 

Harmonic evolutions with incoming gravitational wave: 
Giving the same gravitational wave data (147) as pre¬ 


viously, evolving with the standard boundary condi¬ 
tions (28) and © but using the harmonic Killing gauge 
source functions. It is not obvious how, if at all the space- 
time computed is related to that considered before, but 
in any case we find a very similar decay in ^ 4 . Remark¬ 
ably the growth present in Fig. |10| is absent even in this 
low resolution N = 25 test. 


C. Phasing-in the damped wave gauge 


The transition function: As elsewhere, we follow m 
to transform from one generalized harmonic gauge H\ to 
another H 2 . The composite source function is simply, 

H a (t) = T(t)H 1 a +[l-T(t))H 2 . (152) 

The transition function is, 


T{t) 


0 , t <t t d , 

exp(- {t-t d ) 2 /a 2 d ) , t>t d . 


(153) 
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In the following experiments we choose td = 0 and ad = 
10 M. Note that care must be taken to construct the 
time and space derivatives of H a with the transition func¬ 
tion. This choice results in gauge source functions that 
are only C 1 at t = td, which could be avoided with a 
different transition function. It is not clear if this finite 
differentiability will have a large effect on extracted phys¬ 
ical quantities from a simulation. 

Kerr-Schild, initial slice: For our first phase-in test, 
we started with the Kerr-Schild slicing of Schwarzschild 
and evolved with 74 = 75 = 1/2 and 70 = 1 , on our 
base resolution grid. We took the gauge boundary con¬ 
dition (31) and the constraint preserving condition (39) 


(including a 1/r term). We used the wave gauge param¬ 
eters p = r = 1 and = rjs = 0.1 M. The value of rjg 
here is much smaller than in our wave collapse evolutions. 
The reason for this is that when evolving a blackhole it 
is crucial that the excision boundary is pure outflow in 
the PDEs sense. In other words the characteristic speeds 
must all have the same outward pointing sign. Since the 
speeds in the s l direction are like —/3 s ± a this means 
that the shift can not become too small or else the exci¬ 
sion boundary will fail, which in turn means that r/s can 
not be chosen too large. We therefore place the excision 
boundary deeper into the blackhole so that r m j n = M 
and carefully monitor the coordinate lightspeeds at the 
inner boundary. Note that this requirement is likely to 
cause difficulties when computing extreme gravitational 
waves, because on the one hand large shifts can result in 
poor resolution of important features, but on the other 
they may be required in some other region so that we 
may successfully excise the blackhole region. In the evo¬ 
lution we immediately see significant dynamics and that 
for example the peak of the C x constraint violation along 
the cc-axis is two orders of magnitude greater than in our 
initial Kerr-Schild base run with Killing gauge sources. 
The reason for this is presumably the presence non-trivial 
dynamics, plus the fact that we are excising nearer the 
physical singularity similar to the effect we saw with the 
harmonic Killing slice. Regardless, by t = 100 M the 
data seem very close to stationary. The simulation then 
evolves to the target time t = 1000 M, and remarkably 
at the end of simulation the constraint violation in C x 
along the x-axis has a maximum value which is an order 
of magnitude smaller than in the base run. At no point 
does the excision boundary fail to be outflow. As a check 
of the axisymmetric apparent horizon finder we compare 
the results obtained with the simpler algebraic condition, 


H = 


1 


yj 9rr 


dr log (7 eg) - 2K e g = 0 . 


(154) 


which characterizes the position of the apparent horizon 
in spherical symmetry. We find near perfect agreement 
throughout. The apparent horizon moves from its initial 
radius rjj = 2.00 inwards until it reaches ru = 1-44 
around t = 25. From there the horizon starts to grow 
again and seems to settle down to r# = 1.48. However 
in our lowest resolution run, a small drift of the horizon 



FIG. 12: The radius of the apparent horizon th, and the 
radius at which the outward lightspeed vanishes r c +_ 0 , com¬ 
puted on our base grid with inner boundary at r = 1.2 M. 
To successfully excise, the speed must be negative at the in¬ 
ner boundary. Observe that excision exactly on the apparent 
horizon is not possible throughout all of the run. 


outwards is visible. At late time of the simulation, around 
t = 800, this drift accelerates and we observe that the 
horizon becomes aspherical. Higher resolution runs show 
that this effect converges away. 


Harmonic initial slice: Since the stationary fully har¬ 
monic coordinates are singular at r = M, one might guess 
that the stationary spatial generalized harmonic coordi¬ 
nates with gauge source functions ( 10 ) are also singular 
at some radius on the Killing slice, at least for some range 
of the parameters PLtVs- Given the broad experience in 
using these coordinates in binary blackhole simulations, 
the naive expectation would be that, if present, this co¬ 
ordinate singularity is pushed further towards the physi¬ 
cal singularity rather than out towards the event horizon 
for standard choices of the gauge source functions. But 
this behavior is not clear. To truly resolve the issue one 
could simply solve for such coordinates along the lines 
of [32], but this we defer for the future. Instead we per¬ 
formed simulations varying the initial excision surface 
from the base grid excision radius r m j n = 1.8 M down 
to 7' m in = 1.0 M in steps of 0.2 M. Unsurprisingly we 
find that initially the constraint violation, is greater in 
the excision subpatch as the inner boundary is placed 
closer to the singularity, amounting to about an order of 
magnitude in the C x constraint between the r m j n = M 
and r m j n = 1.2 M boundary runs by t = 50. Besides 
this there is little to distinguish between the five runs, 
and at least down to this excision radius no sign of a 
coordinate singularity forming. By eye the lapse func¬ 
tion in the shared part of the domain agrees very well 
throughout the evolution. Although a slight drift be¬ 
tween them is present towards the end of the test, this is 
acceptable since the outer boundary conditions are being 
imposed at different radii, the solutions need not agree 
everywhere. There is however a time around t = 20 above 
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which the runs with inner boundary r > 1AM fail to be 
outflow at the excision surface. Assuming that this is 
not caused by numerical error this means that bound¬ 
ary conditions are required at the surface. It further¬ 
more means that convergence of the numerical scheme 
as resolution is increased is impossible. The fact that 
this does not correspond to a catastrophic failure of the 
code is inconvenient, because it indicates that great care 
must be taken in monitoring the excision surface. On the 
other hand, since placing the excision boundary very far 
in has a large cost in accuracy, a careful balance must be 
struck. In the SpEC code this is taken care dynamically 
of by a control mechanism lfi.ll 164] which bamps does not 
yet have. In Fig. [12] the relationship between the charac¬ 
ter of the excision boundary and the apparent horizon is 
examined. Comparing the initially harmonic and Kerr- 
Schild slice evolutions with excision radius r m ; n = M we 
find that although the lapse functions initially disagree, 
by about t = 125 M they have exactly the same profile 
and lie almost on top of one another. After this time the 
agreement is maintained. 


VIII. EVOLUTION OF SUPERCRITICAL 
WAVES 


In this section we present the numerical evolution of a 
centered Brill wave, see section V A with A = 8. This 
highly supercritical initial data is used as a test case for 
our excision algorithm for a dynamically forming black- 
hole. 


A. Dynamical excision strategies 


Our dynamical excision method currently consists of 
the following steps: 

1. Evolve to collapse: Evolve on cubed ball grid, 
running the apparent horizon finder in ‘daemon’ mode. 
The finder then triggers a bamps checkpoint once a hori¬ 
zon is found. 

2. Go-to excision grid: Next interpolate the check¬ 
point data onto a cubed-sphere grid. In this step we 
want to place excision boundary just inside the apparent 
horizon, but as we have already seen in the single black- 
hole evolutions this may not always be possible, as some 
wiggle room is needed to allow for dynamical behavior 
of the horizon. This can require some experimentation, 
although fine-tuning does not seem necessary. 

3. Regauge: Adjust the lapse and shift to ensure that 
the excision boundary is pure outflow. As a particular 
choice, we take lapse and shift from Kerr-Schild slicing 
of Schwarzschild, 



( 155 ) 


and translate to Cartesian components in the obvious 
way. It is desirable that the radial coordinate light-speeds 
are close to zero, preferably positive, at the apparent 
horizon, since this determines the direction of motion 
of the horizon. Therefore we choose the to parameter 
to satisfy this condition reasonably well, although again 
without particular fine tuning. 

4. Safety-net evolution: We then use single black- 
hole gauge source parameters like r/^ = 0.1 and rjg = 0.2. 
During the evolution we use a safety net. If any coor¬ 
dinate light-speed on the excision boundary reaches a 
given threshold, typically c* = —0.05 we again regauge 
to guarantee the outflow character is maintained. We 
monitor the apparent horizon, and if it falls off of the 
numerical domain we return to an earlier checkpoint, re¬ 
gauging with a smaller m to avoid this behaviour. As 
the horizon expands we monitor the position and peri¬ 
odically return to the Go-to step above, excising further 
out and regauging with a greater to. 

Discussion: As currently implemented this procedure 
requires that some steps be performed by hand. The 
numerical results in the following subsection serve to 
demonstrate ‘proof of principle’ of this algorithm. On 
the other hand it seems at least clear how those steps 
should be automated. At the regauge step the use of 
the first order GHG variables is again a little fiddly. 
Much more convenient would be if the lapse and shift 
were readily available as variables. But the procedure 
is similar to that described in the gauge perturbation 
tests in section [VII B[ so we do not give full details. Also 
at the regauge step, it might be good to choose lapse 
and shift by abandoning the spherical ansatz and im¬ 
posing that the coordinate light-speeds at the apparent 
horizon vanish. The SpEC approach to controlling the 
excision surface is much more sophisticated, employing 
a control mechanism [63] . we hope to avoid that invest¬ 
ment in the near future. Because we are interested in the 
collapse of waves to form, presumably, a single blackhole 
it seems reasonable to use a simple approach if at all 
possible. One aspect of the method that is not very aes¬ 
thetically appealing, is that by changing the lapse and 
shift in discrete steps we are computing a spacetime, or 
patch of spacetime in coordinates that are not globally 
smooth. Another issue associated with this is that of 
geometric uniqueness, which for the IBVP is an open 
question. Nevertheless one expects that the differences 
to the computed spacetime with one choice of regauging 
parameters or another will be rather small in practice, so 
this does not represent an immediate practical concern. 


B. Supercritical Brill wave evolution 


Initial data and grid setup: We evolved a centered 


Brill wave as described in section V A with seed func¬ 


tion (104). We chose a centered po = 0 wave with A = 8. 
The ADM mass of this initial data is Madm = 1-77. The 
maximum of the Kretschmann scalar in the initial data 
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occurs at the origin, taking the value 1.7 x 10 4 . Follow¬ 
ing the algorithm just outlined, we began on a cubed-ball 
grid with A f cu = 11, AT CS = 13, Af ss = 20, and 55 3 points 
per cube, with internal boundaries r cu = 1.5, r cs = 6.5 
and the outer boundary placed at r = 30 ~ 17 M. 
We ran the code in Cartoon mode on our local clus¬ 
ter Quadler with 240 cores. We evolved with the gen¬ 
eralized harmonic gauge, as in section |VI| in the evolu¬ 
tion of a much weaker A = 2.5 Brill wave, now with the 
gauge parameters vl = 0 and vs = 6. At coordinate 
time t = 1.95 we first found an apparent horizon with 
mass Mh = 1.59 ~ 0.9 M. 

Continuation to code crash: If we continue this evolu¬ 
tion without going to an excision grid after the apparent 
horizon forms, we find that the constraints inside the ap¬ 
parent horizon rapidly grow along with the Kretschmann 
scalar. The run then crashes at roughly t = 3.9. This 
gives the clear signal that if we are to examine the fi¬ 
nal masses of blackholes formed during collapse, using 
the GHG formulation, a robust excision algorithm will 
be essential. In fact at t = 3.85 the horizon has a mass 
of Mh = 1.64 on the cubed-ball grid, but at the end of 
our excision simulation, to be described momentarily, we 
find that 40 M after apparent horizon formation it has 
mass Mh = 1.70. In the first critical gravitational wave 
collapse paper [ 2 ], the blackhole masses were evaluated 
roughly t = 17 M after apparent horizon formation, ac¬ 
cording to a prescription based on the quasinormal modes 
of Schwarzschild. Comparing those values with ours is 
difficult because we use different time coordinates, but 
the basic expectation is that the Maximal slicing condi¬ 
tion is more “singularity avoiding” than one of our gen¬ 
eralized harmonic gauges, and therefore we might expect 
to obtain comparable results if we can evolve for a similar 
coordinate time after the appearance of a horizon. This 
is, however, not clear and deserves further investigation. 
In any case without excising the blackhole region, the 
meager ~ 2M after collapse is clearly insufficient. We 
have seen in i7] that with the moving-puncture method 
this type of data also did not result in successful evolu¬ 
tions beyond apparent horizon formation. But here at 
least a concrete improvement has been made, in that we 
find an apparent horizon before the method fails! 

Evolution on excision grid: Checkpointing the solu¬ 
tion at t = 3.6 we then interpolating, again with barycen- 
tric Lagrange interpolation as used in the apparent hori¬ 
zon finder, onto a cubed-sphere grid with excision radius 
at r = 0.73 M with the outer boundary position fixed, 
and with J\f ss = 27 with 9 angular patches, now with 35 3 
points per cube, naturally again evolving in Cartoon 
mode. In the regauge step we choose here m = 0.4. This 
step immediately removes most of the constraint viola¬ 
tion from the computational domain, and the largest spa¬ 
tial derivatives, so that the constraint monitor is ~ 10~ 8 
as compared to ~ 10 3 on the original cubed-ball. This 
difference seems very troublesome until we take into ac¬ 
count that, for example the peak of the Kretschmann 
scalar on the cubed ball grid is ~ 10 3 whereas on the 


cubed sphere it is ~ 1. So the reduction in the con¬ 
straints obviously occurs because we are removing the 
most extreme part of the domain. Note also that our 
definition of the constraint monitor does not include a 
normalization by the size of the solution, as in for ex¬ 
ample [8j and subsequent papers. In view of this our 
reduction in resolution is justified. The evolution then 
proceeded, now on 120 cores using vl — Vs = 0.1. The 
regauge safety-net was triggered 3 times up to t = 5.9 M, 
having fixed c* = —0.05, but the apparent horizon re¬ 
mains on the computational domain throughout the cal¬ 
culation. At t = 5.9 M we perform the “Go-to” step 
of our algorithm again, this time excising at r = 1.0 M 
choosing m = 0.8. After this the regauge safety-net was 
not called before t = 17 M, when we changed cubed- 
sphere grid once more, keeping the same grid parameters 
but excising at r = 1.12 M, and regauging with m = 1 . 
The evolution continued t = 24.7 M, at which time we 
changed grid for the final time, before which the safety- 
net was again not called. In the last grid we took the exci¬ 
sion radius to be r = 1.24 M and regauged with m = 1.2. 
After this the regauge safety-net was not called, and the 
evolution was terminated at t = 50 M after apparent 
horizon formation. Note that in this evolution the “Go¬ 
to” step also employed the phase-in for the generalized 
harmonic gauge, as described in our single blackhole evo¬ 
lutions in section |VII C[ taking the same parameters em¬ 
ployed in those earlier tests, but now with the initial 
source functions chosen so that the lapse and shift were 
frozen as the evolution starts on the new grid. Other 
experiments show that this procedure is not strictly nec¬ 
essary. It may be that some refinement is required to 
this method to allow the evolution of supercritical data 
indefinitely after the collapse, but examining the mass of 
the apparent horizon, we interpret the solution as having 
mostly settled down, which should be good enough to 
diagnose a final mass of the blackhole. 

Dynamics of the apparent horizon: In the computa¬ 
tion described above, as can be seen in in the left panel 
of Fig. |13[ the apparent horizon is always present on 
the computational domain. The horizon mass initially 
rapidly grows to a value around Mh = 1.7 where it re¬ 
mains roughly constant. Throughout we see that when 
the regauge safety-net is triggered a slight oscillation in 
the horizon mass follows. On the other hand when we 
change grid we see that the horizon mass exhibits a more 
prominent kink. In the right-hand panel of Fig. [13] we 
plot the apparent horizons obtained when, less-wisely, 
the parameter m = 1.4 is chosen in the last “Go-to” 
at t = 24.7 M. With this choice the apparent horizon 
rapidly contracts, although the code fails before it leaves 
the domain. The safety-net is called ever-more frequently 
as the method insists on forcing the inner boundary to 
remain pure outflow, until eventually the code crashes 
at t = 31.6 M. The physical interpretation of this is that 
the excision boundary is falling off of the domain, which 
starts to drift outside the blackhole region, and that the 
safety-net then forces the worldline of the excision bound- 
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FIG. 13: The dynamics of the apparent horizon with our dy¬ 
namical excision strategy for an A = 8 centered Brill wave. 
The green planes indicate the times at which the “Go-to” 
step was applied, and what parameter m was chosen in that 
procedure. The left plot shows a successful choice, and on 
the right what happens if this parameter is chosen less care¬ 
fully. In the upper part of the right hand plot one sees that 
the horizon contracts, and also sees that the ‘regauge’ step is 
frequently applied, resulting in kinks in the horizon. 

ary to be spacelike. This interpretation would be clearer 
if we had an event horizon finder, but is given credence by 
performing evolutions of a Sclrwarzschild blacklrole with 
the m gauge parameter similarly poorly chosen. In such 
tests we see that the areal radius of the excision boundary 
can indeed fall outside of the event horizon at r = 2 M. 


IX. CONCLUSIONS 

We have developed a pseudospectral numerical relativ¬ 
ity code, bamps, and in so doing have made a series of 
improvements and investigations into the approach em¬ 
ployed in the SpEC code. We presented a set of ex¬ 
periments carefully performed so that direct compari¬ 
son with either published work, or independent computa¬ 
tions of the BAM finite differencing code could be made. 
These included evolutions of gauge waves, convergence 
tests, the use of different constraint damping and GHG 
formulation parameters, scaling tests, perturbed black- 
hole evolutions and the treatment of supercritical grav¬ 
itational waves. Ultimately we conclude that the bamps 


code is working efficiently, scales as desired up to large 
numbers of processors, and works on sufficiently general 
grid setups to evolve initial data of interest. Particu¬ 
larly surprising to us was the sensitivity of the method 
to our modifications of the GHG boundary conditions, 
even within the class of constraint preserving conditions. 
This was the case even in our simple evolutions of the 
Schwarzschild spacetime, so it would be very interest¬ 
ing to see the extent to which such results carry over 
to compact binary evolutions, be it in SpEC, or in the 
more distant future in bamps. From the physics point of 
view, however, our focus is presently on the collapse of 
axisymmetric gravitational waves. Much of the develop¬ 
ment reflects this fact. Most notably the implementation 
of octant symmetry with the Cartoon method gives or¬ 
ders of magnitude speedups over evolving the same data 
in full 3d. For a recent complimentary approach see |65j . 
We have additionally written a bespoke axisymmetric ap¬ 
parent horizon finder, which already proved a valuable 
diagnostic tool, crucial in the evolution of supercritical 
data, where the existence of an apparent horizon was 
used as the criterion for moving to an excision grid. 

Naturally further developments to the code may be 
desirable. For physical interpretation, an event horizon 
finder would complement our apparent horizon finder. A 
control system like that of SpEC [63] would be useful in 
controlling the positions of the apparent horizons. But 
the highest priority will likely be in generalizing available 
grid setups to enable dynamical mesh-refinement. 

We have also considered various different types of ax¬ 
isymmetric moment of time-symmetry gravitational wave 
initial data. In forthcoming work we use bamps to evolve 
this initial data, close to the critical amplitude separating 
dispersion and collapse to a blackhole. 
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